Rupture Risk Data Analysis
Ronil V. Chandra
08/07/2020
Objectives:
The purpose of this project is to improve the precision,and understanding of the risk estimate for rupture of a brain aneurysm managed with active surveillance.
This will help clinicians better convey the rupture risks to patients, and is part of the journey toward personalised and patient centred decision making for patients with brain aneurysms.
After performing my systematic review, I have three additional aims:
- To quantitatively synthesise the data in order to calculate a summary pooled proportion.
- To calculate heterogeniety across the studies to help consider whether pooling is appropriate.
- If pooling is appropriate, then to explore the heterogeneity by considering outliers and influential studies, and both pre-specified and post hoc sub-group analysis and meteregression.
The statistical approaches will be tailored to the source dataset and structure of the data which is characterised by the following:
- Rupture of an intracranial aneurysm in a single patient is a rare event, usually estimated at 1% per year or less.
- There are a small number of rupture outcome events in the sample populations, with the proportion very close to zero.
- While the outcome of rupture is binomial, in a sample population, the distribution of rupture outcomes over time is is better described by the non-central hypergeometric distribution, since some aneurysms are considered at higher risk of rupture.
- There are a wide range of size in sample populations identified by the systematic review, with both very small studies and very large studies included.
- The rate of rupture over time does not appear to be constant when biological factors and current published data are considered, thus the risk of rupture will be considered.
To achieve these aims, the various packages and source dataset in R are loaded, and each part carried out with an explanation of the rationale for the data analysis. Version control is carred out with GitHub.
Load required packages
The following R packages were loaded: tidyverse, meta, metafor, BiasedUrn and dmetar.
Load required data
A single excel data file is loaded to carry out all data analysis, with all steps documented below to ensure reproducibility of research.
Correct vector types
Vector types are corrected to integers, factors, doubles as appropriate for analysis in R, and stored as a new tibble.
Calculate individual study proportions and use meta-analytical methods to create a pooled summary proportion
Individual rupture risk at study entry can be calculated by calculating the proportion of patients who ruptured ie pi = xi / ni.
xi = cases ie number of aneurysm ruptures during follow up ni = total ie number of aneurysms at study entry pi = raw proportions ie xi / ni
These can then be combined across studies to consider a meta-analyis of proportions. When considering any meta-analsis, the basic steps are
- calculate a summary statistic for each individual study, such as rupture proportion in our case
- calculate the weighting for each individual study
- choose a random-effect or fixed effect assumption
- calculate the pooled summary statistic
- calculate the p value which communicates the strength of the evidence against the null hypothesis (if there is an intervention group and a control group, but not in our case with a summary proportion)
- derive the confidence intervals which communicates the level of precision or certainty of the summary estimate
- display the meta-analysis result as a forest plot
For every dataset, a suitable effect measure must be chosen, and a choice should be made regarding the meta-analytical methods. Most meta-analytical methods weight the individual effect sizes from each study to create a pooled effect size. In this study, we will consider individual study proportions, and create a pooled summary proportion for rupture risk.
Given unique characteristics of the data that we are synthesising, appropriate statistical methods for meta-analysis and calculation of the confidence intervals must be considered.
Choice of meta-analytical method
The dataset is sparse, and the overall distribution of rupture events across all studies is highly skewed towards zero.
The first step is to improve the statistical properties of a skewed dataset in terms of the data distribution and variance prior to synthesis. Variance is the squared difference from the mean. To achieve this, the data is transformed in order to approximate a normal distribution to enhance the validity of the statistical procedures. The most common transformation methods are the logit and the double arcine transformation.
In the logit transformation, proportions are converted to the natural logarithm of the proportions (i.e., the logit), and assumed to follow the normal distribution. All statistical procedures are performed on the logit proportion, and then the logits converted back to raw proportions for reporting. However, if there are no ruptures and pi = 0, then the logit and variance are undefined. Furthermore, if there is high between-study variablity or small study sample sizes, the logit transformed proportion is underestimated. Hamza 2008 To overcome this statistical limitation, typically a continuity correction of 0.5 is applied. In our analysis, this creates risk of introducing additional sparse data bias and reducing the validity of the result, especially given that pi is close to 0.
The alternate method of transformation to consider which can better normalises the distribution and variance especially for small sample sizes or rare event rates, is the double arcine transformation of Freeman-Tukey Freeman and Tukey 1950. After statistical procedures, the result can be later be backtransformed using the equation derived by Miller Miller 1978. However, the Miller back transformation utilises the harmonic mean of the sample sizes. This affects the backtransformed proportion as described by Schwarzer 2019, and leads to misleading results. This issue is particularly evident when there is a large range of sample sizes. This issue is of particular concern and can lead to misleading results given that the sample sizes on our dataset vary from 22 to 3323.
In the classical meta-analytical methods, after transformation and statistical procedures, the results are backtransformed to study level results, and synthesised.
Data synthesis can be carried out using the generic inverse variance method, which calculates the individual study effect size, and creates a pooled estimate after weighing each study by the inverse of the variance of the effect size estimate. This means that larger studies with smaller standard errors are given more weight than smaller studies with larger standard errors. This minimises the imprecision of the pooled effect size resut. A variation on the generic inverse variance method incorporates an assumption that there is some variation between studies, ie heterogeneity, while measuring the same effect. This produces a random effects meta-analysis, the most common of which is the the DerSimonian and Laird method DerSimonian 1986. This is the most commonly utilised method in medical meta-analysis.
However, when outcome events are rare, such as in our dataset, these classical meta-analytical methods for data-synthesis also have the potential to contribute to additional sparse data bias and give misleading results. Bradburn 2007 from Cochrane
In summary, there are significant statistical limitations of the classical meta-analytical methods of transformation and synthesis with our specific dataset and research question.
To overcome these statistical limitations, we can instead utilise the random intercept logistic regression model for statistical procedures and data synthesis, a type of generalised linear mixed model, as recommended by Stinjen 2010 and Schwarzer 2019.
Our rationale for choice of a GLMM for statistical procedures and data synthesis is based on the following:
- The GLMM is an exact likelihood model based on the the binomial distribution, and thus accounts for the binomial outcome of aneurysm rupture.Hamza 2008
- The GLMM has greater accuracy of the logit proportion and coverage percentages of the corresponding confidence intervals compared to standard normal approximation logit transformation, especially when there are small sample sizes and greater between-study heterogneity. Hamza 2008
- GLMMs can be fitted using the logit transformation, and results backtransformed to the original scale, in our case a proportion.
- The sample populations range from 22 to 3323, which creates misleading results when utilising the Freeman-Tukey methods of transformation and backtransformation of Miller. Schwarzer 2019
- The outcome is very rare, with proportions less than 0.05 is almost all studies, and thus generic inverse variance and DL methods for data synthesis can lead to misleadong results Stinjen 2010
- The GLMM does not require a continuity correction even if pi = 0. Hamza 2008
- Additional calculations to consider covariates can be performed using noncentral hypergeometric models.
The limitation of this statistical approach is that individual study weights utilised to pool the individual studies to create the pooled proportion will not be available.
In the past, utilisation of a GLMM for meta-analysis was not practical due to its computationally intensive nature, and lack of availabilty in standard statistical software packages. However these limitations are now overcome, with statistical modules for GLMMs now available in both SAS and R stastical packages.
Choice of confidence interval
Use of the CI is important, since CIs convey information about magnitude and precision of individual study effect size and the pooled meta-analytical effect size. The choice of the CI should be tailed to the dataset that is present. Options include:
- Wald method or Normal Approximation
- this produces CIs that are often centred on the point estimate
- thus they are not often suitable with proportions close to the boundaries of 0 and 1, since this method can create CIs that are below 0 and above 1.
- for proportions that are often 0 or close to zero, a continuity correction may be applied, but this can lead to additional overshoot
- Clopper-Pearson or Exact method
- Most commonly used, and recommended to avoid approximation by most statistical textbooks.
- Called exact because it is based directly on the binomial distribution and not an approximation.
- Output is usually conservative, and the true coverage is almost always larger than the derived coverage ie closer to a 99% CI than a 95% CI.
- The derived 95% CI does not accurately reflect the true 95% CI unless n is quite large Agresti 2008.
- When n is small eg 10, there is severe overcoverage (closer to 99% coverage) with the true CI much larger than the derived 95% CI.
- Even when n is large, the derived 95% CI does not accurately reflect the true 95% CI when p is near 0 because the binomial distribution is highly skewed. Agresti 1998
- Needs a very large number of n to be accurate Brown 2001
- Wilson method
- Is suggested for small n ie 40 or less and/or extreme probabilities Brown 2001
- The derived CI more accurately reflects the true CI with less variability compared to CP method Agresti 1998
- For small n, the Wilson CI’s are shorter and a more accurate derivation than the CP CI’s Vollset 1993
- Can be used and is recommended for all sample sizes Newcombe 1998
We will choose the Wilson method for CI for the following reasons:
- Rupture of an aneurysm is a rare event with the proportion p very close to zero
- The distribution of the outcomes is highly skewed towards zero, and thus Wald confidence intervals cross zero.
- There is a wide range of sample sizes, with very small studies and very large studies included
- More accurate CIs are likely to be generated using the Wilson method given the highly skewed and sparse dataset.
This is aligned with the recommendations of Vollset 1993, Agresti 1998, Newcombe 1998 and Brown 2001.
Analyse total cohort of all aneurysms 10 mm and less
sizedata <- filter(maindata, size == 10)Assumptions:
Number of patients vs number of aneurysms at entry
Some studies report the number of patients at study entry for a size criteria, while others report the number of aneurysms at study entry, and others report both the number of patients and aneurysms at study entry. Typically the proportion of patients with multiple aneurysms is also reported.
Given that only 1 aneurysm in 1 patient ruptures, which is the outcome, and that the number of aneurysms in that size criteria is most consistently reported, the analysis is performed with the outcome of number of aneurysm ruptures per 100 aneurysms.
This makes biological sense, since aneurysm rupture is considered to be strongly influenced by aneurysm characteristics, 1 patients may have multiple aneurysms with a large range eg 1 to 8 aneurysms, and that conveying a rupture risk per aneurysm is easy to understand from a patient perspective when they harbour multiple aneurysms.
Thus, if the number of aneurysms for a particular size criteria in the study has been reported, this has been extracted and utilised for analysis.
If the number of number of aneuryms is not known for the specific size criteria, and we do know the proportion of patients with multiple aneuryms in the total observation cohort:
- The assumption is that a patient with multiple aneurysms has 2 aneurysms, and we can calculate the total number of aneurysms for that size criteria.
- This can be caculated using the proportion of patients with multiple aneuryms variable multi_tot in the total observation cohort variable num_tot.
- This proportion prop_multi can be multiplied by the number of patients in that specific size criteria. This is the number of patients that we can assume has 2 aneurysms.
- The total number of aneurysms for that size variable total_size is then calculated by adding variable num.
- We have now calculated the total number of aneurysms in that size category, taking into account aneurysm multiplicity.
If the number of number of aneuryms is not known for the specific size criteria, and we do NOT know the proportion of patients with multiple aneuryms in the total observation cohort, but we do know the number of patients in that size category, then 1 patient is assumed to have 1 aneurysm.
For patients with exposure to prior SAH. The number of aneurysms in patients with prior SAH is unknown ie these patients may have multiple aneurysms.
- We can assume that the proportion of patients with multiple aneurysms in the total cohort, is the same as the proportion of patients with multiple aneurysms with prior SAH, and calculate this if we know the number of patients with exposure to prior SAH in the specific size category.
- If the number of patients with exposure to prior SAH in the specific size category is unknown, then we can assume that the proportion of patients with prior SAH in the total observation cohort is the same proportion as those in that size category. Thus if we know the total number of aneurysms in that size category, we can calculate the proportion that have had exposure to prior SAH.
pipe the mutate together *****
dat <- sizedata %>%
mutate(prop_multi = multi_tot / num_tot) %>%
mutate(num_multi = prop_multi * num + num) %>%
mutate(num_multi_temp = coalesce(num_anr, num_multi)) %>%
mutate(total_size_temp = coalesce(num_anr, num)) %>%
mutate(total_size_temp_2 = coalesce(num_multi_temp, total_size_temp)) %>%
mutate(total_size = round(total_size_temp_2, 0)) %>%
mutate(psah_size_temp = psah * prop_multi + psah) %>%
mutate(prop_psah = psah_tot / num_tot) %>%
mutate(num_anr_psah = prop_psah * total_size) %>%
mutate(size_psah_temp = coalesce(psah_size_temp, num_anr_psah)) %>%
mutate(psah_size = round(size_psah_temp, 0)) %>%
unite(auth_year, c(auth, pub), sep = " ", remove = FALSE) %>%
mutate(pop = fct_collapse(sizedata$country, "Japanese" = "Japan", "Non-Japanese" = c("International", "United States", "Switzerland", "Australia", "Korea", "Singapore", "Poland", "China", "Germany", "United Kingdom", "Finland")))Meta-analysis of proportions using GLMM and Wilson CIs with all studies
Now that we have taken into account the number of aneurysms specific to that size category, we can now perform our statistical procedures on that dataset, starting with a meta-analysis of proportions using the GLMM
dat.all <- dat %>%
drop_na(total_size)
pes.summary.glmm.all = metaprop(rupt, total_size,
data=dat.all,
studlab=paste(auth_year),
method="GLMM",
sm="PLOGIT",
method.tau = "ML",
method.ci = "WS",
pscale = 100
)
pes.summary.glmm.all## events 95%-CI
## Bor 2015 0.7444 [ 0.2535; 2.1655]
## Broderick 2009 1.8519 [ 0.5093; 6.5019]
## Burns 2009 0.5780 [ 0.1021; 3.2011]
## Byoun 2016 1.7628 [ 0.9871; 3.1288]
## Choi 2018 0.5780 [ 0.1021; 3.2011]
## Gibbs 2004 0.0000 [ 0.0000; 14.8655]
## Gondar 2016 0.8152 [ 0.2776; 2.3691]
## Guresir 2013 0.7812 [ 0.2660; 2.2715]
## Irazabal 2011 0.0000 [ 0.0000; 7.8652]
## Jeon 2014 0.3521 [ 0.0966; 1.2746]
## Jiang 2013 0.0000 [ 0.0000; 7.1348]
## Juvela 2013 18.6747 [13.4796; 25.2868]
## Loumiotis 2011 0.0000 [ 0.0000; 2.3446]
## Matsubara 2004 0.0000 [ 0.0000; 2.3736]
## Matsumoto 2013 2.6549 [ 0.9069; 7.5160]
## Mizoi 1995 0.0000 [ 0.0000; 15.4639]
## Morita 2012 1.8658 [ 1.4582; 2.3845]
## Murayama 2016 2.2384 [ 1.6660; 3.0014]
## Oh 2013 0.0000 [ 0.0000; 16.8179]
## Serrone 2016 0.5155 [ 0.0911; 2.8615]
## So 2010 1.1450 [ 0.3902; 3.3118]
## Sonobe 2010 1.5625 [ 0.7589; 3.1897]
## Teo 2016 2.3810 [ 0.8130; 6.7666]
## Thien 2017 0.6173 [ 0.1090; 3.4133]
## Tsukahara 2005 3.4722 [ 1.4921; 7.8703]
## Tsutsumi 2000 4.3478 [ 1.4896; 12.0212]
## Villablanca 2013 1.5544 [ 0.5300; 4.4697]
## Wiebers-R 1998 1.3503 [ 0.8558; 2.1244]
## Wilkinson 2018 0.0000 [ 0.0000; 14.8655]
## Zylkowski 2015 2.7273 [ 0.9318; 7.7131]
##
## Number of studies combined: k = 30
##
## events 95%-CI
## Fixed effect model 1.7872 [1.5638; 2.0419]
## Random effects model 1.2112 [0.7844; 1.8659]
##
## Quantifying heterogeneity:
## tau^2 = 0.7974; tau = 0.8930; I^2 = 82.7%; H = 2.41
##
## Test of heterogeneity:
## Q d.f. p-value Test
## 170.38 29 < 0.0001 Wald-type
## 148.41 29 < 0.0001 Likelihood-Ratio
##
## Details on meta-analytical method:
## - Random intercept logistic regression model
## - Maximum-likelihood estimator for tau^2
## - Logit transformation
## - Wilson Score confidence interval for individual studies
## - Continuity correction of 0.5 in studies with zero cell frequencies
## (only used to calculate individual study results)
## - Events per 100 observationsThe output from the random effects meta-analysis using the GLMM:
Summary estimate is 1.2112 [0.7844; 1.8659]
Note that for GLMMs no continuity correction is used. Meta documentation Given our rationale for choosing the GLMM, this should produce the least biased result and reasonable coverage probabilities for the 95% CI, as suggested by Stinjen 2010. Note CIs are using Wilson score method.
Remember that the confidence interval from this random-effects meta-analysis describes uncertainty in the location of average proportion across the individual studies. Thus there is likely to be an even higher degree of uncertainty in the true population rupture risk. Cochrane handbook 10.10.4.2
Display meta-analysis result using a publication quality forest plot
The easiest way to communicate the result of the statistical procedure is via a forest plot.
forest(pes.summary.glmm.all,
layout = "meta",
comb.fixed = FALSE,
comb.random = TRUE,
leftlabs = c("Study", "Ruptures", "Total"),
rightcols = c("effect", "ci"),
rightlabs = c("Ruptures per 100", "95% CI"),
smlab = " ",
xlim=c(0,10),
xlab = "Rupture Proportion per 100",
pooled.events = TRUE,
) 
Heterogeniety
What is heterogeniety
Heterogeniety refers to all types of variation across studies. Heterogeniety can be considered in terms of
- Clinical heterogeneity, which refers to variablity in characteristics of patients, exposures, interventions and outcomes.
- Methodological heterogeneity, which refers to variablity in study designs eg randomization, blinding etc.
- Statistical heterogeneity, which refers to the situation when variation across studies is greater than variation or random error within a study.
Clinical and methodolocal heterogeniety cannot be calculated. These are both subjectively evaluated by the meta-analyst. If clinical and methodological differences are not significant, then quantative synthesis is considered appropriate.
If there is substantial clinical and/or methdological heterogeneity, a quantitative synthesis should not be performed to create a pooled effect measure. A qualitative or narrative systematic review should be performed. This is because, the observed effects across studies are greater than what is expected due to random chance alone, which limits the generalisability of the result. The analogy is that a systematic review brings together apples and oranges, and that combining these in the setting of high heterogeneity yields a meaningless result.
Since statistical heterogeneity is a consequence of clinical and/or methodological heterogeniety, this always occurs in a meta-analysis, statistical heterogeneity is inevitable Higgins et al 2003
How to identify and measure statistical heterogeniety
Overall, identification of heterogeniety is critical to assist in interpreting the results of the meta-analyis. This gives the reader confidence that the effect demonstrated is generalisable. It is impossible to completely avoid heterogeneity, however clinical and methodological heterogeneity can be minimised by using strict inclusion criteria during systematic review.
A simple way to identify heterogeneity is to review the forest plot, and see if the confidence intervals overlap. If there is poor overlap, there is statistical heterogeneity.
Statitical heterogeneity can also be measured using Cochrane’s Q statistic, tau-squared statistic or the Inconsistency measure I^2 statistic.
- Cochrane’s Q. This test performs poorly for small numbers of studies. If the p value for the Q statistic is below 0.05, then this suggests that there is significant between-study heterogeneity.
- tau-squared statistic. Uncommonly reported. An estimate of the between-study variance. When tau-squared is zero this is indicative of no heterogeneity.
- I^2 Inconsistency statistic. I^2 represents the proportion of total heterogeniety that can be attributed to the actual between-study heterogeniety, and can be classified into the degree of heterogeneity.
From Cochrane:
0% to 40%: might not be important 30% to 60%: may represent moderate heterogeneity 50% to 90%: may represent substantial heterogeneity 75% to 100%: considerable heterogeneity
95% CIs can also be calculated for the I^2 statistic to demonstrate the uncertainty of the result.
We will choose the I^2 statistic for the following reasons:
- This measure is not affected by the number of studies included.
- The output offers a percentage that is easily interpreted
- This is the standard method of reporting for Cochrane reviews.
pes.summary.glmm.all## events 95%-CI
## Bor 2015 0.7444 [ 0.2535; 2.1655]
## Broderick 2009 1.8519 [ 0.5093; 6.5019]
## Burns 2009 0.5780 [ 0.1021; 3.2011]
## Byoun 2016 1.7628 [ 0.9871; 3.1288]
## Choi 2018 0.5780 [ 0.1021; 3.2011]
## Gibbs 2004 0.0000 [ 0.0000; 14.8655]
## Gondar 2016 0.8152 [ 0.2776; 2.3691]
## Guresir 2013 0.7812 [ 0.2660; 2.2715]
## Irazabal 2011 0.0000 [ 0.0000; 7.8652]
## Jeon 2014 0.3521 [ 0.0966; 1.2746]
## Jiang 2013 0.0000 [ 0.0000; 7.1348]
## Juvela 2013 18.6747 [13.4796; 25.2868]
## Loumiotis 2011 0.0000 [ 0.0000; 2.3446]
## Matsubara 2004 0.0000 [ 0.0000; 2.3736]
## Matsumoto 2013 2.6549 [ 0.9069; 7.5160]
## Mizoi 1995 0.0000 [ 0.0000; 15.4639]
## Morita 2012 1.8658 [ 1.4582; 2.3845]
## Murayama 2016 2.2384 [ 1.6660; 3.0014]
## Oh 2013 0.0000 [ 0.0000; 16.8179]
## Serrone 2016 0.5155 [ 0.0911; 2.8615]
## So 2010 1.1450 [ 0.3902; 3.3118]
## Sonobe 2010 1.5625 [ 0.7589; 3.1897]
## Teo 2016 2.3810 [ 0.8130; 6.7666]
## Thien 2017 0.6173 [ 0.1090; 3.4133]
## Tsukahara 2005 3.4722 [ 1.4921; 7.8703]
## Tsutsumi 2000 4.3478 [ 1.4896; 12.0212]
## Villablanca 2013 1.5544 [ 0.5300; 4.4697]
## Wiebers-R 1998 1.3503 [ 0.8558; 2.1244]
## Wilkinson 2018 0.0000 [ 0.0000; 14.8655]
## Zylkowski 2015 2.7273 [ 0.9318; 7.7131]
##
## Number of studies combined: k = 30
##
## events 95%-CI
## Fixed effect model 1.7872 [1.5638; 2.0419]
## Random effects model 1.2112 [0.7844; 1.8659]
##
## Quantifying heterogeneity:
## tau^2 = 0.7974; tau = 0.8930; I^2 = 82.7%; H = 2.41
##
## Test of heterogeneity:
## Q d.f. p-value Test
## 170.38 29 < 0.0001 Wald-type
## 148.41 29 < 0.0001 Likelihood-Ratio
##
## Details on meta-analytical method:
## - Random intercept logistic regression model
## - Maximum-likelihood estimator for tau^2
## - Logit transformation
## - Wilson Score confidence interval for individual studies
## - Continuity correction of 0.5 in studies with zero cell frequencies
## (only used to calculate individual study results)
## - Events per 100 observationsI^2 = 82.7%, which is high, which means that the between study variation due to clinical and methodological variation between studies is high, and greater than that expected by random chance variation within studies.
This may be explained by co-variates and thus consideration should be made to avoid data-sythesis or expore reasons for heterogeneity.
For interest, the tau-squared = 0.9098 and the Q-statistic is < 0.0001 which also confirms high between-study heterogeniety.
We could have identified the heterogeneity by reviewing the overlap between confidence intervals in the forest plot, as we can see that there is a clear outlier - Juvela et al.
How to explore and address heterogeneity
There are various strategies to explore and adddress heterogeneity.
Check for any manual data entry errors. If data has been entered into the statistical software incorrectly, this can result in inaccurate results. By utilising a single excel file as the source of data, with no additional manual data entry, this risk is minimised.
Check that a random-effects meta-analysis was performed. Clinical and methodological heterogeneity is always present in medical studies, and in our case a random-effects model has ben used.
Avoid quantitative data synthesis and pooling studies. If clinical and methodological heterogeneity are substantial, then a quantitative data synthesis shuld not be peformed. A narrative or qualitative review should be performed instead.
Choosing to pool the data, while exploring heterogeneity. Heterogeneity can be explored through sub-group analysis or metaregression to help identify study-level effect modifiers. This is where the outcome varies in different population or intervention characteristics. Such sub-group and metaregression analyses are best be pre-specified, and regardless they shoud be interpreted with caution and considered hypothesis generating.
Exclude studies. 1 or 2 studies may be outliers or have high influence on the pooled effect size. In general, exclusion of studies should be avoided, since it may introduce bias. However, if there is an obvious reason for the outlying result, then the study can be excluded with confidence. Overall, it is often best to perform the analysis both with and without the outlier / influential study as part of a sensitivity analysis.
Investigating heterogeneity
The forest plot reveals an outlier with the Juvela et al study. To commence our investigation into heterogeneity, we will first consider outlier and influence analysis.
Outlier studies with extreme effect sizes can cause increased between study heterogenieity. This may arise from small or low quality studies. In addition, overly influential studies can push the pooled effect size higher or lower, which means that the pooled result is less robust since it relies on 1 or a small number of studies.
There are a number of ways to detect outliers and influential studies.
Outlier assessment
One way to identify statistical outliers is to see if the confidence intervals of one study do not overlap with the pooled effect size. This outlier has an extreme effect size, which means that it cannot be included in the total study population pooled to create the pooled effect size. This is because inclusion of the statistical outlier reduces the robustness of the pooled effect size result.
A common way to detect an outlier are to identify studies where
- The upper bound of the 95% CI is lower than the lower bound of the pooled effect 95% CI (i.e., extremely small effects)
- The lower bound of the 95% CI is higher than the upper bound of the pooled effect CI (i.e., extremely large effects)
The current results reveal that that the pooled rupture risk is 1.2112 [0.7844; 1.8659]. So lets identify studies that are outliers.
find.outliers(pes.summary.glmm.all)## Identified outliers (fixed-effect model)
## ----------------------------------------
## "Jeon 2014", "Juvela 2013"
##
## Results with outliers removed
## -----------------------------
## events 95%-CI exclude
## Bor 2015 0.7444 [ 0.2535; 2.1655]
## Broderick 2009 1.8519 [ 0.5093; 6.5019]
## Burns 2009 0.5780 [ 0.1021; 3.2011]
## Byoun 2016 1.7628 [ 0.9871; 3.1288]
## Choi 2018 0.5780 [ 0.1021; 3.2011]
## Gibbs 2004 0.0000 [ 0.0000; 14.8655]
## Gondar 2016 0.8152 [ 0.2776; 2.3691]
## Guresir 2013 0.7812 [ 0.2660; 2.2715]
## Irazabal 2011 0.0000 [ 0.0000; 7.8652]
## Jeon 2014 0.3521 [ 0.0966; 1.2746] *
## Jiang 2013 0.0000 [ 0.0000; 7.1348]
## Juvela 2013 18.6747 [13.4796; 25.2868] *
## Loumiotis 2011 0.0000 [ 0.0000; 2.3446]
## Matsubara 2004 0.0000 [ 0.0000; 2.3736]
## Matsumoto 2013 2.6549 [ 0.9069; 7.5160]
## Mizoi 1995 0.0000 [ 0.0000; 15.4639]
## Morita 2012 1.8658 [ 1.4582; 2.3845]
## Murayama 2016 2.2384 [ 1.6660; 3.0014]
## Oh 2013 0.0000 [ 0.0000; 16.8179]
## Serrone 2016 0.5155 [ 0.0911; 2.8615]
## So 2010 1.1450 [ 0.3902; 3.3118]
## Sonobe 2010 1.5625 [ 0.7589; 3.1897]
## Teo 2016 2.3810 [ 0.8130; 6.7666]
## Thien 2017 0.6173 [ 0.1090; 3.4133]
## Tsukahara 2005 3.4722 [ 1.4921; 7.8703]
## Tsutsumi 2000 4.3478 [ 1.4896; 12.0212]
## Villablanca 2013 1.5544 [ 0.5300; 4.4697]
## Wiebers-R 1998 1.3503 [ 0.8558; 2.1244]
## Wilkinson 2018 0.0000 [ 0.0000; 14.8655]
## Zylkowski 2015 2.7273 [ 0.9318; 7.7131]
##
## Number of studies combined: k = 28
##
## events 95%-CI
## Fixed effect model 1.6086 [1.3908; 1.8598]
## Random effects model 1.4042 [1.0703; 1.8401]
##
## Quantifying heterogeneity:
## tau^2 = 0.0680; tau = 0.2608; I^2 = 26.9%; H = 1.17
##
## Test of heterogeneity:
## Q d.f. p-value Test
## 22.09 27 0.7326 Wald-type
## 41.85 27 0.0340 Likelihood-Ratio
##
## Details on meta-analytical method:
## - Random intercept logistic regression model
## - Maximum-likelihood estimator for tau^2
## - Logit transformation
## - Wilson Score confidence interval for individual studies
## - Continuity correction of 0.5 in studies with zero cell frequencies
## (only used to calculate individual study results)
## - Events per 100 observations
##
## Identified outliers (random-effects model)
## ------------------------------------------
## "Juvela 2013"
##
## Results with outliers removed
## -----------------------------
## events 95%-CI exclude
## Bor 2015 0.7444 [ 0.2535; 2.1655]
## Broderick 2009 1.8519 [ 0.5093; 6.5019]
## Burns 2009 0.5780 [ 0.1021; 3.2011]
## Byoun 2016 1.7628 [ 0.9871; 3.1288]
## Choi 2018 0.5780 [ 0.1021; 3.2011]
## Gibbs 2004 0.0000 [ 0.0000; 14.8655]
## Gondar 2016 0.8152 [ 0.2776; 2.3691]
## Guresir 2013 0.7812 [ 0.2660; 2.2715]
## Irazabal 2011 0.0000 [ 0.0000; 7.8652]
## Jeon 2014 0.3521 [ 0.0966; 1.2746]
## Jiang 2013 0.0000 [ 0.0000; 7.1348]
## Juvela 2013 18.6747 [13.4796; 25.2868] *
## Loumiotis 2011 0.0000 [ 0.0000; 2.3446]
## Matsubara 2004 0.0000 [ 0.0000; 2.3736]
## Matsumoto 2013 2.6549 [ 0.9069; 7.5160]
## Mizoi 1995 0.0000 [ 0.0000; 15.4639]
## Morita 2012 1.8658 [ 1.4582; 2.3845]
## Murayama 2016 2.2384 [ 1.6660; 3.0014]
## Oh 2013 0.0000 [ 0.0000; 16.8179]
## Serrone 2016 0.5155 [ 0.0911; 2.8615]
## So 2010 1.1450 [ 0.3902; 3.3118]
## Sonobe 2010 1.5625 [ 0.7589; 3.1897]
## Teo 2016 2.3810 [ 0.8130; 6.7666]
## Thien 2017 0.6173 [ 0.1090; 3.4133]
## Tsukahara 2005 3.4722 [ 1.4921; 7.8703]
## Tsutsumi 2000 4.3478 [ 1.4896; 12.0212]
## Villablanca 2013 1.5544 [ 0.5300; 4.4697]
## Wiebers-R 1998 1.3503 [ 0.8558; 2.1244]
## Wilkinson 2018 0.0000 [ 0.0000; 14.8655]
## Zylkowski 2015 2.7273 [ 0.9318; 7.7131]
##
## Number of studies combined: k = 29
##
## events 95%-CI
## Fixed effect model 1.5475 [1.3390; 1.7879]
## Random effects model 1.2625 [0.9445; 1.6858]
##
## Quantifying heterogeneity:
## tau^2 = 0.1359; tau = 0.3687; I^2 = 41.8%; H = 1.31
##
## Test of heterogeneity:
## Q d.f. p-value Test
## 27.45 28 0.4940 Wald-type
## 49.86 28 0.0067 Likelihood-Ratio
##
## Details on meta-analytical method:
## - Random intercept logistic regression model
## - Maximum-likelihood estimator for tau^2
## - Logit transformation
## - Wilson Score confidence interval for individual studies
## - Continuity correction of 0.5 in studies with zero cell frequencies
## (only used to calculate individual study results)
## - Events per 100 observationsHere using the random effects model, Juvela is identified as the outlier.
The pooled rupture risk changes with Juvela et al is 1.2112 [0.7844; 1.8659]
The pooled rupture risk without Juvela et al is 1.2625 [0.9445; 1.6858].
Although the pooled rupture risk is still between 1-2%, the precision of the estimate has improved, and the I2 has reduced from 82 to 42%, which makes the result more robust and more generalisable.
Influence analysis
Apart from excluding statistical outliers which reduce the robustness of the pooled result, it is important to identify studies in the meta-analysis that extert high influence on the pooled results. If 1 study very highly influences the pooled effect size result, this also makes the pooled effect size result less robust and less generalisable.
A common method to identify influential studies is the leave-1-out method, where the results of the meta-analysis are recalculated leaving out 1 study at a time. This helps identify individual studies that may exert very high influence on the result, pushing the pooled effect size higher or lower.
inf.analysis <- InfluenceAnalysis(x = pes.summary.glmm.all,
random = TRUE,
subplot.heights = c(30,18),
subplot.widths = c(30,30),
forest.lims = c(0,0.03))## [===========================================================================] DONEsummary(inf.analysis)## Leave-One-Out Analysis (Sorted by I2)
## -----------------------------------
## Effect LLCI ULCI I2
## Omitting Juvela 2013 -4.359 -4.653 -4.066 0.418
## Omitting Morita 2012 -4.441 -4.905 -3.976 0.793
## Omitting Murayama 2016 -4.448 -4.911 -3.986 0.803
## Omitting Jeon 2014 -4.336 -4.773 -3.899 0.821
## Omitting Loumiotis 2011 -4.342 -4.771 -3.914 0.823
## Omitting Matsubara 2004 -4.343 -4.772 -3.914 0.823
## Omitting Wiebers-R 1998 -4.421 -4.886 -3.956 0.827
## Omitting Jiang 2013 -4.376 -4.813 -3.939 0.830
## Omitting Irazabal 2011 -4.378 -4.815 -3.941 0.831
## Omitting Serrone 2016 -4.370 -4.814 -3.927 0.831
## Omitting Burns 2009 -4.375 -4.820 -3.931 0.832
## Omitting Byoun 2016 -4.433 -4.897 -3.970 0.832
## Omitting Choi 2018 -4.375 -4.820 -3.931 0.832
## Omitting Gibbs 2004 -4.388 -4.827 -3.950 0.832
## Omitting Mizoi 1995 -4.389 -4.827 -3.950 0.832
## Omitting Oh 2013 -4.390 -4.829 -3.951 0.832
## Omitting Thien 2017 -4.378 -4.823 -3.933 0.832
## Omitting Tsukahara 2005 -4.452 -4.907 -3.996 0.832
## Omitting Tsutsumi 2000 -4.449 -4.902 -3.997 0.832
## Omitting Wilkinson 2018 -4.388 -4.827 -3.950 0.832
## Omitting Bor 2015 -4.381 -4.833 -3.928 0.833
## Omitting Guresir 2013 -4.383 -4.837 -3.930 0.833
## Omitting Gondar 2016 -4.386 -4.840 -3.932 0.834
## Omitting Matsumoto 2013 -4.437 -4.894 -3.980 0.836
## Omitting So 2010 -4.404 -4.861 -3.947 0.836
## Omitting Sonobe 2010 -4.425 -4.887 -3.962 0.836
## Omitting Zylkowski 2015 -4.438 -4.894 -3.981 0.836
## Omitting Broderick 2009 -4.420 -4.875 -3.965 0.837
## Omitting Teo 2016 -4.433 -4.891 -3.976 0.837
## Omitting Villablanca 2013 -4.418 -4.876 -3.959 0.837
##
##
## Influence Diagnostics
## -------------------
## rstudent dffits cook.d cov.r QE.del hat weight
## Omitting Bor 2015 -0.855 -0.173 0.030 1.040 171.145 0.039 3.919
## Omitting Broderick 2009 0.117 0.042 0.002 1.061 175.212 0.033 3.299
## Omitting Burns 2009 -0.849 -0.131 0.017 1.022 173.290 0.023 2.272
## Omitting Byoun 2016 0.089 0.053 0.003 1.099 174.272 0.053 5.320
## Omitting Choi 2018 -0.849 -0.131 0.017 1.022 173.290 0.023 2.272
## Omitting Gibbs 2004 0.175 0.030 0.001 1.025 175.338 0.014 1.375
## Omitting Gondar 2016 -0.754 -0.148 0.022 1.046 171.774 0.039 3.918
## Omitting Guresir 2013 -0.801 -0.160 0.026 1.043 171.485 0.039 3.919
## Omitting Irazabal 2011 -0.262 -0.025 0.001 1.022 175.028 0.014 1.387
## Omitting Jeon 2014 -1.553 -0.325 0.102 0.987 167.836 0.033 3.322
## Omitting Jiang 2013 -0.328 -0.033 0.001 1.022 174.939 0.014 1.388
## Omitting Juvela 2013 10.287 -0.420 0.016 0.218 30.914 0.058 5.769
## Omitting Loumiotis 2011 -1.054 -0.131 0.017 1.007 173.227 0.014 1.396
## Omitting Matsubara 2004 -1.046 -0.129 0.017 1.007 173.253 0.014 1.396
## Omitting Matsumoto 2013 0.512 0.124 0.016 1.068 175.305 0.039 3.891
## Omitting Mizoi 1995 0.203 0.033 0.001 1.025 175.341 0.014 1.374
## Omitting Morita 2012 0.170 0.077 0.006 1.111 170.016 0.060 5.990
## Omitting Murayama 2016 0.406 0.132 0.019 1.106 175.135 0.059 5.918
## Omitting Oh 2013 0.264 0.040 0.002 1.025 175.341 0.014 1.372
## Omitting Serrone 2016 -0.943 -0.148 0.022 1.018 172.945 0.023 2.273
## Omitting So 2010 -0.385 -0.061 0.004 1.063 173.676 0.039 3.913
## Omitting Sonobe 2010 -0.058 0.015 0.000 1.090 174.030 0.049 4.939
## Omitting Teo 2016 0.395 0.102 0.011 1.071 175.342 0.039 3.895
## Omitting Thien 2017 -0.796 -0.121 0.015 1.024 173.475 0.023 2.271
## Omitting Tsukahara 2005 0.870 0.203 0.042 1.066 174.594 0.046 4.553
## Omitting Tsutsumi 2000 1.048 0.216 0.047 1.048 174.213 0.039 3.866
## Omitting Villablanca 2013 -0.058 0.011 0.000 1.071 174.777 0.039 3.907
## Omitting Wiebers-R 1998 -0.247 -0.031 0.001 1.098 168.872 0.056 5.619
## Omitting Wilkinson 2018 0.175 0.030 0.001 1.025 175.338 0.014 1.375
## Omitting Zylkowski 2015 0.541 0.129 0.017 1.068 175.284 0.039 3.890
## infl
## Omitting Bor 2015
## Omitting Broderick 2009
## Omitting Burns 2009
## Omitting Byoun 2016
## Omitting Choi 2018
## Omitting Gibbs 2004
## Omitting Gondar 2016
## Omitting Guresir 2013
## Omitting Irazabal 2011
## Omitting Jeon 2014
## Omitting Jiang 2013
## Omitting Juvela 2013
## Omitting Loumiotis 2011
## Omitting Matsubara 2004
## Omitting Matsumoto 2013
## Omitting Mizoi 1995
## Omitting Morita 2012
## Omitting Murayama 2016
## Omitting Oh 2013
## Omitting Serrone 2016
## Omitting So 2010
## Omitting Sonobe 2010
## Omitting Teo 2016
## Omitting Thien 2017
## Omitting Tsukahara 2005
## Omitting Tsutsumi 2000
## Omitting Villablanca 2013
## Omitting Wiebers-R 1998
## Omitting Wilkinson 2018
## Omitting Zylkowski 2015
##
##
## Baujat Diagnostics (sorted by Heterogeneity Contribution)
## -------------------------------------------------------
## HetContrib InfluenceEffectSize
## Omitting Juvela 2013 126.815 34.314
## Omitting Jeon 2014 7.434 13.094
## Omitting Wiebers-R 1998 5.914 13.006
## Omitting Bor 2015 4.137 15.146
## Omitting Guresir 2013 3.801 15.336
## Omitting Morita 2012 3.759 14.157
## Omitting Gondar 2016 3.517 15.496
## Omitting Serrone 2016 2.386 16.290
## Omitting Loumiotis 2011 2.110 16.151
## Omitting Matsubara 2004 2.084 16.171
## Omitting Burns 2009 2.043 16.504
## Omitting Choi 2018 2.043 16.504
## Omitting Thien 2017 1.858 16.617
## Omitting So 2010 1.642 16.571
## Omitting Sonobe 2010 1.268 16.625
## Omitting Tsutsumi 2000 1.114 18.594
## Omitting Byoun 2016 1.014 16.786
## Omitting Tsukahara 2005 0.731 18.822
## Omitting Villablanca 2013 0.557 17.284
## Omitting Jiang 2013 0.402 17.278
## Omitting Irazabal 2011 0.314 17.330
## Omitting Murayama 2016 0.165 19.577
## Omitting Broderick 2009 0.129 17.672
## Omitting Zylkowski 2015 0.057 18.157
## Omitting Matsumoto 2013 0.037 18.125
## Omitting Gibbs 2004 0.004 17.569
## Omitting Wilkinson 2018 0.004 17.569
## Omitting Mizoi 1995 0.001 17.580
## Omitting Oh 2013 0.001 17.600
## Omitting Teo 2016 0.000 17.988plot(inf.analysis, "baujat")
The Baujat Plot (Baujat et al. 2002) is a diagnostic plot to detect studies overly contributing to the heterogeneity. The study contribution to Cochran’s Q is plotted on the horizontal axis, and the study influence on the pooled effect size on the vertical axis. All studies on the right side of the plot are of interest, since the contribute to heterogeniety.
Here Juvela et al is identified as both a high contributor to heterogeneity, and high influence on the pooled result.
plot(inf.analysis, "influence")
We can also perform additional plots to assess influence of each study. These plots proposed by Viechtbauer & Cheung (2010) also demonstrate that the Juvela study is not only a statistical outlier, but also an influential study. This is important because this further corroborates that the Juvela study adds statistical heterogeneity, and influences the overall pooled result.
plot(inf.analysis, "i2")
In these 2 plots, we can see the impact of the leave 1 out analysis on the pooled effect size. This is ordered by the reduction in statistical heterogeneity as measured by I2.
We can again see that the precision of the pooled effect size estimate improves with exlcusion of Juvela et al, and that the I2 reduces the most with exclusion of Juvela et al.
These findings corroborate the findings of the outlier analysis, the Baujat plots, and influence analysis proposed by Viechtbauer and Cheung that Juvela is the main source of heterogeniety.
Why Juvela et al excluded.
The outlier and influence analysis are concordant that Juvela et al is an outlier which influences the pooled effect estimate and reduces the precision of the pooled effect estimate. This is on the basis of statistical heterogeneity, and that the rupture proportion observed in the Juvela et al study varies greater than what is expected due to random chance alone.
The analogy is that while the remaining studies are different varieties of apples, that Juvela study is an orange, and that combining this study with all other studies will yield much less meaningful result.
Since statistical heterogeneity is a consequence of clinical and/or methodological heterogeniety, we can explore if there are particular clinical or methodogical factors that are likely to be responsible.
Proposed factors to explore are proportion of patients with exposure to prior subarachnoid haemorrhage, patient enrollment period, and length of follow up.
Meta-analysis of proportions using GLMM and Wilson CIs with Juvela excluded
Here is the main study result, with the Juvela study excluded due to reasons previously stated.
dat <- dat %>%
slice(-12) %>%
drop_na(total_size)
pes.summary.glmm = metaprop(rupt, total_size,
data=dat,
studlab=paste(auth_year),
method="GLMM",
sm="PLOGIT",
method.tau = "ML",
method.ci = "WS",
pscale = 100
)
pes.summary.glmm## events 95%-CI
## Bor 2015 0.7444 [0.2535; 2.1655]
## Broderick 2009 1.8519 [0.5093; 6.5019]
## Burns 2009 0.5780 [0.1021; 3.2011]
## Byoun 2016 1.7628 [0.9871; 3.1288]
## Choi 2018 0.5780 [0.1021; 3.2011]
## Gibbs 2004 0.0000 [0.0000; 14.8655]
## Gondar 2016 0.8152 [0.2776; 2.3691]
## Guresir 2013 0.7812 [0.2660; 2.2715]
## Irazabal 2011 0.0000 [0.0000; 7.8652]
## Jeon 2014 0.3521 [0.0966; 1.2746]
## Jiang 2013 0.0000 [0.0000; 7.1348]
## Loumiotis 2011 0.0000 [0.0000; 2.3446]
## Matsubara 2004 0.0000 [0.0000; 2.3736]
## Matsumoto 2013 2.6549 [0.9069; 7.5160]
## Mizoi 1995 0.0000 [0.0000; 15.4639]
## Morita 2012 1.8658 [1.4582; 2.3845]
## Murayama 2016 2.2384 [1.6660; 3.0014]
## Oh 2013 0.0000 [0.0000; 16.8179]
## Serrone 2016 0.5155 [0.0911; 2.8615]
## So 2010 1.1450 [0.3902; 3.3118]
## Sonobe 2010 1.5625 [0.7589; 3.1897]
## Teo 2016 2.3810 [0.8130; 6.7666]
## Thien 2017 0.6173 [0.1090; 3.4133]
## Tsukahara 2005 3.4722 [1.4921; 7.8703]
## Tsutsumi 2000 4.3478 [1.4896; 12.0212]
## Villablanca 2013 1.5544 [0.5300; 4.4697]
## Wiebers-R 1998 1.3503 [0.8558; 2.1244]
## Wilkinson 2018 0.0000 [0.0000; 14.8655]
## Zylkowski 2015 2.7273 [0.9318; 7.7131]
##
## Number of studies combined: k = 29
##
## events 95%-CI
## Fixed effect model 1.5475 [1.3390; 1.7879]
## Random effects model 1.2625 [0.9445; 1.6858]
##
## Quantifying heterogeneity:
## tau^2 = 0.1359; tau = 0.3687; I^2 = 41.8%; H = 1.31
##
## Test of heterogeneity:
## Q d.f. p-value Test
## 27.45 28 0.4940 Wald-type
## 49.86 28 0.0067 Likelihood-Ratio
##
## Details on meta-analytical method:
## - Random intercept logistic regression model
## - Maximum-likelihood estimator for tau^2
## - Logit transformation
## - Wilson Score confidence interval for individual studies
## - Continuity correction of 0.5 in studies with zero cell frequencies
## (only used to calculate individual study results)
## - Events per 100 observationsThe output from the random effects meta-analysis using the GLMM:
Summary estimate is 1.2625 [0.9445; 1.6858] with I2 of 41.8%.
Display meta-analysis result using a publication quality forest plot
The easiest way to communicate the result of the statistical procedure is via a forest plot.
forest(pes.summary.glmm,
layout = "meta",
comb.fixed = FALSE,
comb.random = TRUE,
leftlabs = c("Study", "Ruptures", "Total"),
rightcols = c("effect", "ci"),
rightlabs = c("Ruptures per 100", "95% CI"),
smlab = " ",
xlim=c(0,10),
xlab = "Rupture Proportion per 100",
pooled.events = TRUE,
) 
Exploring residual heterogeneity.
Now that the pooled effect size is more homogenous, we can explore residual heterogeniety through subgroup analysis and meta-regression.
Subgroup analyses require splitting of the data into subgroups to make comparison. Typically this data is hard to extract from the published literature, and thus are better utilised when individual patient level data are available for pooling. Nonetheless, we will perform the sub-group meta-analysis above to explore residual heterogeneity and consider the limitations of the analyses.
Meta-regression allows the effects of multiple continuous and categorical variables to be investigated simultaneously, unlike subgroup analysis which considers one categorical variable only. This generally utilises a random effect model to conduct the analysis in each subgroup, and then considers additional statistical testing to compare the pooled results across the subgroups. Meta-regression should generally not be considered when there are less than ten studies in a meta-analysis. The meta-regression cooeffient obtained describes how the outcome (dependent variable) changes withn unit increase in the potential effect modifier. If the outcome is a ratio measure, the log-transformed value should be used in the regression model. Cochrane Handbook 10.11.4
Meta-regression and sub-group analysis are non-randomized and observational in nature. It may not be possible to ajust for all potential confounders and thus they can can lead to misleading conclusions. Even if there are a large number of covariates adjusted for, we cannot be certain that all potential confounders have been identified. These analyses cannot prove causality, and at best, they can be considered hypothesis generating.
This is an important concept, because if subgroup analysis or metaregression findings are presented as definitive conclusions there is risk of people being denied an effective intervention or treated with an ineffective or harmful intervention. This can also generate misleading recommendations about future research directions, which if followed would waste scarce research resources.
Moving forward, we shall choose to pool the data, and explore heterogeneity in the synthesised data with sub-group analyses and meta-regression.
Subgroup analysis
The following sub-group group analyses will be performed.
Pre-specified sub-group analysis of rupture risk by size of aneurysm (7mm and less, 5mm and less, and 3mm and less).
Prespecified sub-group analysis of rupture risk by size of aneurysm by exposure to prior subarachnoid haemorrhage.
or study proportion with PSAH?
Post hoc sub-group analysis of rupture risk comparing prospective studies to retrospective studies.
Post hoc sub-group analysis of rupture risk comparing study source population from Japan vs non-Japanese popultions.
Post hoc sub-group analysis of rupture risk with length of follow up categorised as intermediate (5 years or less), or long term (>5 years)
Subgroup analysis - aneurysm size categories
sizedata7 <- filter(maindata, size == 7) dat7 <- sizedata7 %>%
mutate(prop_multi = multi_tot / num_tot) %>%
mutate(num_multi = prop_multi * num + num) %>%
mutate(num_multi_temp = coalesce(num_anr, num_multi)) %>%
mutate(total_size_temp = coalesce(num_anr, num)) %>%
mutate(total_size_temp_2 = coalesce(num_multi_temp, total_size_temp)) %>%
mutate(total_size = round(total_size_temp_2, 0)) %>%
mutate(psah_size_temp = psah * prop_multi + psah) %>%
mutate(prop_psah = psah_tot / num_tot) %>%
mutate(num_anr_psah = prop_psah * total_size) %>%
mutate(size_psah_temp = coalesce(psah_size_temp, num_anr_psah)) %>%
mutate(psah_size = round(size_psah_temp, 0)) %>%
unite(auth_year, c(auth, pub), sep = " ", remove = FALSE) %>%
mutate(pop = fct_collapse(sizedata$country,
"Japanese" = "Japan",
"Non-Japanese" = c("International", "United States", "Switzerland", "Australia", "Korea", "Singapore", "Poland", "China", "Germany", "United Kingdom", "Finland")
))Meta-analysis of proportions of 7 mm aneurysms and less using GLMM and Wilson CIs
dat7 <- dat7 %>%
slice(-12) %>%
drop_na(total_size)
pes.summary.glmm7 = metaprop(rupt, total_size,
data=dat7,
studlab=paste(auth_year),
method="GLMM",
sm="PLOGIT",
method.tau = "ML",
method.ci = "WS",
pscale = 100
)
pes.summary.glmm7## events 95%-CI
## Bor 2015 0.4963 [0.1362; 1.7912]
## Broderick 2009 1.8519 [0.5093; 6.5019]
## Byoun 2016 1.8613 [1.0424; 3.3018]
## Choi 2018 0.5780 [0.1021; 3.2011]
## Gibbs 2004 0.0000 [0.0000; 14.8655]
## Gondar 2016 0.8152 [0.2776; 2.3691]
## Guresir 2013 0.7812 [0.2660; 2.2715]
## Irazabal 2011 0.0000 [0.0000; 8.2010]
## Jeon 2014 0.3521 [0.0966; 1.2746]
## Jiang 2013 0.0000 [0.0000; 8.7622]
## Matsubara 2004 0.0000 [0.0000; 2.9815]
## Matsumoto 2013 0.8850 [0.1564; 4.8431]
## Mizoi 1995 0.0000 [0.0000; 15.4639]
## Morita 2012 1.2933 [0.9397; 1.7774]
## Murayama 2016 2.0679 [1.5164; 2.8142]
## Oh 2013 0.0000 [0.0000; 21.5311]
## Serrone 2016 0.0000 [0.0000; 1.9417]
## So 2010 0.4695 [0.0829; 2.6110]
## Sonobe 2010 1.5625 [0.7589; 3.1897]
## Teo 2016 2.3810 [0.8130; 6.7666]
## Thien 2017 0.0000 [0.0000; 2.5637]
## Tsukahara 2005 1.8692 [0.5141; 6.5604]
## Tsutsumi 2000 3.3333 [0.9189; 11.3638]
## Villablanca 2013 1.5544 [0.5300; 4.4697]
## Wiebers-P 2003 0.8580 [0.4520; 1.6225]
## Wilkinson 2018 0.0000 [0.0000; 16.1125]
## Zylkowski 2015 2.7273 [0.9318; 7.7131]
##
## Number of studies combined: k = 27
##
## events 95%-CI
## Fixed effect model 1.2613 [1.0638; 1.4950]
## Random effects model 1.0580 [0.7710; 1.4503]
##
## Quantifying heterogeneity:
## tau^2 = 0.1416; tau = 0.3763; I^2 = 37.5%; H = 1.26
##
## Test of heterogeneity:
## Q d.f. p-value Test
## 23.23 26 0.6201 Wald-type
## 43.60 26 0.0167 Likelihood-Ratio
##
## Details on meta-analytical method:
## - Random intercept logistic regression model
## - Maximum-likelihood estimator for tau^2
## - Logit transformation
## - Wilson Score confidence interval for individual studies
## - Continuity correction of 0.5 in studies with zero cell frequencies
## (only used to calculate individual study results)
## - Events per 100 observationsPublication quality forest plots for 7mm and less
forest(pes.summary.glmm7,
layout = "meta",
comb.fixed = FALSE,
comb.random = TRUE,
leftlabs = c("Study", "Ruptures", "Total"),
rightcols = c("effect", "ci"),
rightlabs = c("Ruptures per 100", "95% CI"),
smlab = " ",
xlim=c(0,10),
xlab = "Rupture Proportion per 100",
pooled.events = TRUE,
) 
sub-group analysis based on aneurysm size 5 mm and less
sizedata5 <- filter(maindata, size == 5) dat5 <- sizedata5 %>%
mutate(prop_multi = multi_tot / num_tot) %>%
mutate(num_multi = prop_multi * num + num) %>%
mutate(num_multi_temp = coalesce(num_anr, num_multi)) %>%
mutate(total_size_temp = coalesce(num_anr, num)) %>%
mutate(total_size_temp_2 = coalesce(num_multi_temp, total_size_temp)) %>%
mutate(total_size = round(total_size_temp_2, 0)) %>%
mutate(psah_size_temp = psah * prop_multi + psah) %>%
mutate(prop_psah = psah_tot / num_tot) %>%
mutate(num_anr_psah = prop_psah * total_size) %>%
mutate(size_psah_temp = coalesce(psah_size_temp, num_anr_psah)) %>%
mutate(psah_size = round(size_psah_temp, 0)) %>%
unite(auth_year, c(auth, pub), sep = " ", remove = FALSE) %>%
mutate(pop = fct_collapse(sizedata$country,
"Japanese" = "Japan",
"Non-Japanese" = c("International", "United States", "Switzerland", "Australia", "Korea", "Singapore", "Poland", "China", "Germany", "United Kingdom", "Finland")
))Subgroup meta-analysis of proportions of 5 mm aneurysms and less using GLMM and Wilson CIs
dat5 <- dat5 %>%
slice(-12) %>%
drop_na(total_size) %>%
drop_na(rupt)
pes.summary.glmm5 = metaprop(rupt, total_size,
data=dat5,
studlab=paste(auth_year),
method="GLMM",
sm="PLOGIT",
method.tau = "ML",
method.ci = "WS",
pscale = 100
)
pes.summary.glmm5## events 95%-CI
## Bor 2015 0.7435 [0.2041; 2.6699]
## Broderick 2009 1.3889 [0.2456; 7.4566]
## Gibbs 2004 0.0000 [0.0000; 16.8179]
## Gondar 2016 0.9091 [0.3096; 2.6383]
## Guresir 2013 1.0563 [0.3599; 3.0592]
## Irazabal 2011 0.0000 [0.0000; 10.1515]
## Jeon 2014 0.3521 [0.0966; 1.2746]
## Jiang 2013 0.0000 [0.0000; 16.8179]
## Matsubara 2004 0.0000 [0.0000; 2.9815]
## Matsumoto 2013 1.2658 [0.2238; 6.8276]
## Mizoi 1995 0.0000 [0.0000; 15.4639]
## Morita 2012 1.1500 [0.7675; 1.7198]
## Murayama 2016 1.2813 [0.8477; 1.9325]
## Oh 2013 0.0000 [0.0000; 22.8095]
## Serrone 2016 0.0000 [0.0000; 7.4100]
## So 2010 0.4695 [0.0829; 2.6110]
## Sonobe 2010 1.5625 [0.7589; 3.1897]
## Thien 2017 0.0000 [0.0000; 2.5637]
## Tsukahara 2005 1.8692 [0.5141; 6.5604]
## Tsutsumi 2000 3.3333 [0.9189; 11.3638]
## Wilkinson 2018 0.0000 [0.0000; 16.8179]
## Zylkowski 2015 2.1739 [0.5982; 7.5835]
##
## Number of studies combined: k = 22
##
## events 95%-CI
## Fixed effect model 1.0624 [0.8427; 1.3385]
## Random effects model 1.0624 [0.8427; 1.3385]
##
## Quantifying heterogeneity:
## tau^2 = 0; tau = 0; I^2 = 0.0%; H = 1.00
##
## Test of heterogeneity:
## Q d.f. p-value Test
## 8.53 21 0.9924 Wald-type
## 19.59 21 0.5475 Likelihood-Ratio
##
## Details on meta-analytical method:
## - Random intercept logistic regression model
## - Maximum-likelihood estimator for tau^2
## - Logit transformation
## - Wilson Score confidence interval for individual studies
## - Continuity correction of 0.5 in studies with zero cell frequencies
## (only used to calculate individual study results)
## - Events per 100 observationsPublication quality forest plots for 5mm and less
forest(pes.summary.glmm5,
layout = "meta",
comb.fixed = FALSE,
comb.random = TRUE,
leftlabs = c("Study", "Ruptures", "Total"),
rightcols = c("effect", "ci"),
rightlabs = c("Ruptures per 100", "95% CI"),
smlab = " ",
xlim=c(0,10),
xlab = "Rupture Proportion per 100",
pooled.events = TRUE,
) 
sub-group analysis based on aneurysm size 3 mm and less
sizedata3 <- filter(maindata, size == 3) dat3 <- sizedata3 %>%
mutate(prop_multi = multi_tot / num_tot) %>%
mutate(num_multi = prop_multi * num + num) %>%
mutate(num_multi_temp = coalesce(num_anr, num_multi)) %>%
mutate(total_size_temp = coalesce(num_anr, num)) %>%
mutate(total_size_temp_2 = coalesce(num_multi_temp, total_size_temp)) %>%
mutate(total_size = round(total_size_temp_2, 0)) %>%
mutate(psah_size_temp = psah * prop_multi + psah) %>%
mutate(prop_psah = psah_tot / num_tot) %>%
mutate(num_anr_psah = prop_psah * total_size) %>%
mutate(size_psah_temp = coalesce(psah_size_temp, num_anr_psah)) %>%
mutate(psah_size = round(size_psah_temp, 0)) %>%
unite(auth_year, c(auth, pub), sep = " ", remove = FALSE) %>%
mutate(pop = fct_collapse(sizedata$country,
"Japanese" = "Japan",
"Non-Japanese" = c("International", "United States", "Switzerland", "Australia", "Korea", "Singapore", "Poland", "China", "Germany", "United Kingdom", "Finland")
))subgroup meta-analysis of proportions of 3 mm aneurysms and less using GLMM and Wilson CIs
dat3 <- dat3 %>%
slice(-12) %>%
drop_na(total_size) %>%
drop_na(rupt)
pes.summary.glmm3 = metaprop(rupt, total_size,
data=dat3,
studlab=paste(auth_year),
method="GLMM",
control=list(stepadj=0.5,maxiter=1000),
sm="PLOGIT",
method.tau = "ML",
method.ci = "WS",
pscale = 100
)
pes.summary.glmm3## events 95%-CI
## Bor 2015 3.0303 [0.5369; 15.3187]
## Broderick 2009 1.3889 [0.2456; 7.4566]
## Gibbs 2004 0.0000 [0.0000; 24.2494]
## Gondar 2016 0.5376 [0.0950; 2.9821]
## Guresir 2013 1.6000 [0.4399; 5.6463]
## Irazabal 2011 0.0000 [0.0000; 15.4639]
## Jiang 2013 0.0000 [0.0000; 16.8179]
## Matsumoto 2013 0.0000 [0.0000; 4.6371]
## Mizoi 1995 0.0000 [0.0000; 15.4639]
## Oh 2013 0.0000 [0.0000; 56.1497]
## Serrone 2016 0.0000 [0.0000; 7.4100]
## So 2010 1.6393 [0.2900; 8.7189]
## Sonobe 2010 1.6129 [0.4434; 5.6903]
## Wilkinson 2018 0.0000 [0.0000; 17.5879]
## Zylkowski 2015 0.0000 [0.0000; 11.6970]
##
## Number of studies combined: k = 15
##
## events 95%-CI
## Fixed effect model 0.9401 [0.4708; 1.8683]
## Random effects model 0.9401 [0.4708; 1.8683]
##
## Quantifying heterogeneity:
## tau^2 = 0; tau = 0; I^2 = 0.0%; H = 1.00
##
## Test of heterogeneity:
## Q d.f. p-value Test
## 1.60 14 1.0000 Wald-type
## 7.46 14 0.9156 Likelihood-Ratio
##
## Details on meta-analytical method:
## - Random intercept logistic regression model
## - Maximum-likelihood estimator for tau^2
## - Logit transformation
## - Wilson Score confidence interval for individual studies
## - Continuity correction of 0.5 in studies with zero cell frequencies
## (only used to calculate individual study results)
## - Events per 100 observationsPublication quality forest plots for 3mm and less
forest(pes.summary.glmm3,
layout = "meta",
comb.fixed = FALSE,
comb.random = TRUE,
leftlabs = c("Study", "Ruptures", "Total"),
rightcols = c("effect", "ci"),
rightlabs = c("Ruptures per 100", "95% CI"),
smlab = " ",
xlim=c(0,10),
xlab = "Rupture Proportion per 100",
pooled.events = TRUE,
) 
Sub-group analysis based on study type - retrospective compared to prospective
pes.summary.glmm.type = metaprop(rupt, total_size,
data=dat,
studlab=paste(auth_year),
byvar = type,
bylab = "Study Type",
method="GLMM",
sm="PLOGIT",
method.tau = "ML",
method.ci = "WS",
pscale = 100
)
pes.summary.glmm.type## events 95%-CI Study Type
## Bor 2015 0.7444 [0.2535; 2.1655] prospective
## Broderick 2009 1.8519 [0.5093; 6.5019] prospective
## Burns 2009 0.5780 [0.1021; 3.2011] retrospective
## Byoun 2016 1.7628 [0.9871; 3.1288] retrospective
## Choi 2018 0.5780 [0.1021; 3.2011] retrospective
## Gibbs 2004 0.0000 [0.0000; 14.8655] retrospective
## Gondar 2016 0.8152 [0.2776; 2.3691] prospective
## Guresir 2013 0.7812 [0.2660; 2.2715] retrospective
## Irazabal 2011 0.0000 [0.0000; 7.8652] retrospective
## Jeon 2014 0.3521 [0.0966; 1.2746] retrospective
## Jiang 2013 0.0000 [0.0000; 7.1348] retrospective
## Loumiotis 2011 0.0000 [0.0000; 2.3446] prospective
## Matsubara 2004 0.0000 [0.0000; 2.3736] retrospective
## Matsumoto 2013 2.6549 [0.9069; 7.5160] retrospective
## Mizoi 1995 0.0000 [0.0000; 15.4639] retrospective
## Morita 2012 1.8658 [1.4582; 2.3845] prospective
## Murayama 2016 2.2384 [1.6660; 3.0014] prospective
## Oh 2013 0.0000 [0.0000; 16.8179] retrospective
## Serrone 2016 0.5155 [0.0911; 2.8615] retrospective
## So 2010 1.1450 [0.3902; 3.3118] retrospective
## Sonobe 2010 1.5625 [0.7589; 3.1897] prospective
## Teo 2016 2.3810 [0.8130; 6.7666] retrospective
## Thien 2017 0.6173 [0.1090; 3.4133] retrospective
## Tsukahara 2005 3.4722 [1.4921; 7.8703] retrospective
## Tsutsumi 2000 4.3478 [1.4896; 12.0212] retrospective
## Villablanca 2013 1.5544 [0.5300; 4.4697] retrospective
## Wiebers-R 1998 1.3503 [0.8558; 2.1244] retrospective
## Wilkinson 2018 0.0000 [0.0000; 14.8655] retrospective
## Zylkowski 2015 2.7273 [0.9318; 7.7131] retrospective
##
## Number of studies combined: k = 29
##
## events 95%-CI
## Fixed effect model 1.5475 [1.3390; 1.7879]
## Random effects model 1.2625 [0.9445; 1.6858]
##
## Quantifying heterogeneity:
## tau^2 = 0.1359; tau = 0.3687; I^2 = 41.8%; H = 1.31
##
## Quantifying residual heterogeneity:
## I^2 = 0.0% [0.0%; 39.7%]; H = 1.00 [1.00; 1.29]
##
## Test of heterogeneity:
## Q d.f. p-value Test
## 27.45 28 0.4940 Wald-type
## 49.86 28 0.0067 Likelihood-Ratio
##
## Results for subgroups (fixed effect model):
## k events 95%-CI Q I^2
## Study Type = prospective 7 1.7828 [1.4927; 2.1280] 6.33 38.6%
## Study Type = retrospective 22 1.2286 [0.9571; 1.5759] 19.68 24.4%
##
## Test for subgroup differences (fixed effect model):
## Q d.f. p-value
## Between groups 5.69 1 0.0170
## Within groups 26.01 27 0.5178
##
## Results for subgroups (random effects model):
## k events 95%-CI tau^2 tau
## Study Type = prospective 7 1.5673 [0.9451; 2.5885] 0.0530 0.2302
## Study Type = retrospective 22 1.1668 [0.8213; 1.6552] 0.1230 0.3507
##
## Test for subgroup differences (random effects model):
## Q d.f. p-value
## Between groups 0.89 1 0.3466
##
## Details on meta-analytical method:
## - Random intercept logistic regression model
## - Maximum-likelihood estimator for tau^2
## - Logit transformation
## - Wilson Score confidence interval for individual studies
## - Continuity correction of 0.5 in studies with zero cell frequencies
## (only used to calculate individual study results)
## - Events per 100 observationsThere is greater imprecision in the result for prospective studies with higher heterogeniety. Retrospective studies have less heterogeneity and narrower confidence intervals, but are subject to greater degrees of bias.
forest(pes.summary.glmm.type,
layout = "meta",
comb.fixed = FALSE,
comb.random = TRUE,
leftlabs = c("Study", "Ruptures", "Total"),
rightcols = c("effect", "ci"),
rightlabs = c("Ruptures per 100", "95% CI"),
smlab = " ",
xlim=c(0,10),
xlab = "Rupture Proportion per 100",
pooled.events = TRUE,
) 
Sub-group analysis based on study population - Japanese vs non-Japanese
pes.summary.glmm.pop = metaprop(rupt, total_size,
data=dat,
studlab=paste(auth_year),
byvar = pop,
bylab = "Population",
method="GLMM",
sm="PLOGIT",
method.tau = "ML",
method.ci = "WS",
pscale = 100
)
pes.summary.glmm.pop## events 95%-CI Population
## Bor 2015 0.7444 [0.2535; 2.1655] Non-Japanese
## Broderick 2009 1.8519 [0.5093; 6.5019] Non-Japanese
## Burns 2009 0.5780 [0.1021; 3.2011] Non-Japanese
## Byoun 2016 1.7628 [0.9871; 3.1288] Non-Japanese
## Choi 2018 0.5780 [0.1021; 3.2011] Non-Japanese
## Gibbs 2004 0.0000 [0.0000; 14.8655] Non-Japanese
## Gondar 2016 0.8152 [0.2776; 2.3691] Non-Japanese
## Guresir 2013 0.7812 [0.2660; 2.2715] Non-Japanese
## Irazabal 2011 0.0000 [0.0000; 7.8652] Non-Japanese
## Jeon 2014 0.3521 [0.0966; 1.2746] Non-Japanese
## Jiang 2013 0.0000 [0.0000; 7.1348] Non-Japanese
## Loumiotis 2011 0.0000 [0.0000; 2.3446] Non-Japanese
## Matsubara 2004 0.0000 [0.0000; 2.3736] Japanese
## Matsumoto 2013 2.6549 [0.9069; 7.5160] Japanese
## Mizoi 1995 0.0000 [0.0000; 15.4639] Japanese
## Morita 2012 1.8658 [1.4582; 2.3845] Japanese
## Murayama 2016 2.2384 [1.6660; 3.0014] Japanese
## Oh 2013 0.0000 [0.0000; 16.8179] Non-Japanese
## Serrone 2016 0.5155 [0.0911; 2.8615] Non-Japanese
## So 2010 1.1450 [0.3902; 3.3118] Non-Japanese
## Sonobe 2010 1.5625 [0.7589; 3.1897] Japanese
## Teo 2016 2.3810 [0.8130; 6.7666] Non-Japanese
## Thien 2017 0.6173 [0.1090; 3.4133] Non-Japanese
## Tsukahara 2005 3.4722 [1.4921; 7.8703] Non-Japanese
## Tsutsumi 2000 4.3478 [1.4896; 12.0212] Japanese
## Villablanca 2013 1.5544 [0.5300; 4.4697] Non-Japanese
## Wiebers-R 1998 1.3503 [0.8558; 2.1244] Non-Japanese
## Wilkinson 2018 0.0000 [0.0000; 14.8655] Non-Japanese
## Zylkowski 2015 2.7273 [0.9318; 7.7131] Non-Japanese
##
## Number of studies combined: k = 29
##
## events 95%-CI
## Fixed effect model 1.5475 [1.3390; 1.7879]
## Random effects model 1.2625 [0.9445; 1.6858]
##
## Quantifying heterogeneity:
## tau^2 = 0.1359; tau = 0.3687; I^2 = 41.8%; H = 1.31
##
## Quantifying residual heterogeneity:
## I^2 = 0.0% [0.0%; 25.0%]; H = 1.00 [1.00; 1.15]
##
## Test of heterogeneity:
## Q d.f. p-value Test
## 27.45 28 0.4940 Wald-type
## 49.86 28 0.0067 Likelihood-Ratio
##
## Results for subgroups (fixed effect model):
## k events 95%-CI Q I^2
## Population = Non-Japanese 22 1.1164 [0.8731; 1.4266] 17.54 14.6%
## Population = Japanese 7 1.9494 [1.6300; 2.3300] 3.35 0.0%
##
## Test for subgroup differences (fixed effect model):
## Q d.f. p-value
## Between groups 12.97 1 0.0003
## Within groups 20.89 27 0.7912
##
## Results for subgroups (random effects model):
## k events 95%-CI tau^2 tau
## Population = Non-Japanese 22 1.0665 [0.7710; 1.4735] 0.0634 0.2518
## Population = Japanese 7 1.9494 [1.6300; 2.3300] 0 0
##
## Test for subgroup differences (random effects model):
## Q d.f. p-value
## Between groups 10.25 1 0.0014
##
## Details on meta-analytical method:
## - Random intercept logistic regression model
## - Maximum-likelihood estimator for tau^2
## - Logit transformation
## - Wilson Score confidence interval for individual studies
## - Continuity correction of 0.5 in studies with zero cell frequencies
## (only used to calculate individual study results)
## - Events per 100 observationsforest(pes.summary.glmm.pop,
layout = "meta",
comb.fixed = FALSE,
comb.random = TRUE,
leftlabs = c("Study", "Ruptures", "Total"),
rightcols = c("effect", "ci"),
rightlabs = c("Ruptures per 100", "95% CI"),
smlab = " ",
xlim=c(0,10),
xlab = "Rupture Proportion per 100",
pooled.events = TRUE,
) 
Sub-group meta-analysis of proportions for risk of rupture in patients with and without exposure to prior SAH
how about changing categorising this as proportion of patients in study with prior SAH at <10% vs 10% or more? 9 studies vs other - yes ****************
Firstly structure the data in individual studies with authors in a 2 x 2 table
ai = aneuryms with exposure to prior SAH and rupture +ve = variable rupt_psah bi = aneuryms with exposure to prior SAH and rupture -ve ci = aneuryms with no exposure to prior SAH, and rupture +ve = rupt minus rupt_psah di = aneuryms with no exposure to prior SAH, and rupture -ve
n1i = ai + bi = total aneurysms with prior SAH n2i = ci + di = total aneurysms without prior SAH
n1i + n2i = total_anr included with and without prior SAH
total number of events = total number of ruptures = variable rupt
Assumptions:
- Only 1 aneurysm ruptures per patient.
- If number of patients with prior SAH is known, then this is multiplied by the multiplicity proportion to calculate number of aneurysms with exposure to prior SAH in the specific size criteria.
- If the patients of patients with prior is SAH for the size criteria is unknown, then the proportion of patients with prior SAH in the total observation cohort of all aneurysms is carried forward and applied to the number of aneurysms, and an estimated number of aneurysms with prior SAH is calculated based on the total observation cohort.
- Data analysis is limited to studies in which there are known or calculated values of ai, bi, ci, and di using the assumptions. If there are NAs in the 2 x 2 table, then that study is excluded.
dat.psah <- dat %>%
unite(auth_year, c(auth, pub), sep = " ", remove = FALSE) %>%
drop_na(total_size) %>%
rename(ai = rupt_psah) %>%
rename(rupt.total = rupt) %>%
mutate(ci = rupt.total - ai) %>%
rename(n1i = psah_size) %>%
mutate(n2i = total_size - n1i) %>%
mutate(bi = n1i - ai) %>%
mutate(di = n2i - ci) %>%
select(auth_year, ai, bi, ci, di, n1i, n2i) %>%
drop_na(ai, bi, ci, di) %>%
mutate_if(is.numeric, round, 0)then create new tibble for patients with and without history of prior subarachnoid haemorrhage.
dat.psahpos <- dat.psah %>%
filter(n1i!=0)
dat.psahneg <- dat.psah %>%
filter(n2i!=0) run new GLMM (random intercept logistic regression model) for patients with prior history of subaarachnoid haemorrhage. Studies that did not include any patients with prior history of subarachnoid haemorrage are excluded.
pes.summary.psahpos = metaprop(ai, n1i,
data=dat.psahpos,
studlab=paste(auth_year),
method="GLMM",
sm="PLOGIT",
method.tau = "ML",
method.ci = "WS",
pscale = 100
)
pes.summary.psahpos## events 95%-CI
## Gondar 2016 4.7619 [0.8456; 22.6694]
## Matsubara 2004 0.0000 [0.0000; 39.0334]
## Mizoi 1995 0.0000 [0.0000; 65.7620]
## Serrone 2016 0.0000 [0.0000; 65.7620]
## Sonobe 2010 2.0833 [0.3687; 10.8992]
## Teo 2016 0.0000 [0.0000; 8.7622]
## Tsukahara 2005 0.0000 [0.0000; 11.3513]
## Wiebers-R 1998 2.1197 [1.3276; 3.3683]
## Wilkinson 2018 0.0000 [0.0000; 56.1497]
## Zylkowski 2015 2.8571 [0.5061; 14.5331]
##
## Number of studies combined: k = 10
##
## events 95%-CI
## Fixed effect model 2.0222 [1.3083; 3.1136]
## Random effects model 2.0222 [1.3083; 3.1136]
##
## Quantifying heterogeneity:
## tau^2 = 0; tau = 0; I^2 = 0.0%; H = 1.00
##
## Test of heterogeneity:
## Q d.f. p-value Test
## 0.70 9 0.9999 Wald-type
## 4.12 9 0.9035 Likelihood-Ratio
##
## Details on meta-analytical method:
## - Random intercept logistic regression model
## - Maximum-likelihood estimator for tau^2
## - Logit transformation
## - Wilson Score confidence interval for individual studies
## - Continuity correction of 0.5 in studies with zero cell frequencies
## (only used to calculate individual study results)
## - Events per 100 observationsforest(pes.summary.psahpos,
layout = "meta",
comb.fixed = FALSE,
comb.random = TRUE,
leftlabs = c("Study", "Ruptures", "Total"),
rightcols = c("effect", "ci"),
rightlabs = c("Ruptures per 100", "95% CI"),
smlab = " ",
xlim=c(0,30),
xlab = "Rupture Proportion per 100",
pooled.events = TRUE,
) 
Then run new GLMM (random intercept logistic regression model) for patients without history of PSAH
pes.summary.psahneg = metaprop(ci, n2i,
data=dat.psahneg,
studlab=paste(auth_year),
method="GLMM",
sm="PLOGIT",
method.tau = "ML",
method.ci = "WS",
pscale = 100
)
pes.summary.psahneg## events 95%-CI
## Broderick 2009 1.8519 [0.5093; 6.5019]
## Gibbs 2004 0.0000 [0.0000; 14.8655]
## Gondar 2016 0.5764 [0.1582; 2.0768]
## Guresir 2013 0.7812 [0.2660; 2.2715]
## Irazabal 2011 0.0000 [0.0000; 7.8652]
## Jeon 2014 0.3521 [0.0966; 1.2746]
## Jiang 2013 0.0000 [0.0000; 7.1348]
## Loumiotis 2011 0.0000 [0.0000; 2.3446]
## Matsubara 2004 0.0000 [0.0000; 2.4650]
## Mizoi 1995 0.0000 [0.0000; 16.8179]
## Oh 2013 0.0000 [0.0000; 16.8179]
## Serrone 2016 0.5208 [0.0920; 2.8907]
## Sonobe 2010 1.5000 [0.6892; 3.2335]
## Teo 2016 3.4884 [1.1934; 9.7609]
## Thien 2017 0.6173 [0.1090; 3.4133]
## Tsukahara 2005 4.3860 [1.8878; 9.8581]
## Tsutsumi 2000 4.3478 [1.4896; 12.0212]
## Wiebers-R 1998 0.1883 [0.0333; 1.0589]
## Wilkinson 2018 0.0000 [0.0000; 16.8179]
## Zylkowski 2015 2.6667 [0.7344; 9.2115]
##
## Number of studies combined: k = 20
##
## events 95%-CI
## Fixed effect model 0.8802 [0.6197; 1.2489]
## Random effects model 0.7882 [0.4301; 1.4404]
##
## Quantifying heterogeneity:
## tau^2 = 0.6924; tau = 0.8321; I^2 = 53.2%; H = 1.46
##
## Test of heterogeneity:
## Q d.f. p-value Test
## 25.77 19 0.1368 Wald-type
## 37.22 19 0.0074 Likelihood-Ratio
##
## Details on meta-analytical method:
## - Random intercept logistic regression model
## - Maximum-likelihood estimator for tau^2
## - Logit transformation
## - Wilson Score confidence interval for individual studies
## - Continuity correction of 0.5 in studies with zero cell frequencies
## (only used to calculate individual study results)
## - Events per 100 observationsHere we can see that the I^2 = 53.2% within the sub-group, which means that there remains moderate between-study variation, even within the subgroup. In other words, even if the study falls within this subgroup, we cannot accurately estimate the pooled rupture risk.
A forest plot is constructed for rupture risk in patients without history of prior SAH below.
forest(pes.summary.psahneg,
layout = "meta",
comb.fixed = FALSE,
comb.random = TRUE,
leftlabs = c("Study", "Ruptures", "Total"),
rightcols = c("effect", "ci"),
rightlabs = c("Ruptures per 100", "95% CI"),
smlab = " ",
xlim=c(0,10),
xlab = "Rupture Proportion per 100",
pooled.events = TRUE,
) 
Subgroup analyis for length of follow up
A clinical hypothesis to explore is that the proportion of patients that rupture may be associated with the length of follow up in the study.
This exploratory analysis considers the proportion of ruptures that occur in studies with intermediate length of follow up defined as 5 years or less, and long term follow up defined as more than 5 years. 5 Years is chosen as it is a relatively common measure for patient level discussions about risk vs benefit across a number of fields of medicine.
dat.fu <- dat %>%
mutate(fu = coalesce(fu_mean_tot,fu_med_tot)) %>%
drop_na(total_size) %>%
unite(auth_year, c(auth, pub), sep = " ", remove = FALSE) %>%
select(auth_year, rupt, total_size, fu) Calculate mean and median follow up in included studies
dat.fu.mean <- dat.fu %>%
drop_na(fu) %>%
summarise(mean = mean(fu, na.rm = TRUE))
dat.fu.mean## # A tibble: 1 x 1
## mean
## <dbl>
## 1 43.3dat.fu.median <- dat.fu %>%
drop_na(fu) %>%
summarise(median = median(fu, na.rm = TRUE))
dat.fu.median## # A tibble: 1 x 1
## median
## <dbl>
## 1 38.4Mean follow up length across included studies is 43 months, and median is 38 months. It seems more appropriate to utilise the means as a summary compared to the median.
Now that the column with length of follow up is created, we can perform a subgroup meta-analysis of studies with intermediate and long term follow up.
should we consider short term follow up vs longer term follow up categorised at 3 years given that median was 3 years? ***********
dat.fu.total <- dat.fu %>%
mutate(fu_5yr = case_when(fu <= 60 ~ "Less than 5 Years", fu > 60 ~ "Greater than 5 years")) %>%
drop_na(fu_5yr)pes.summary.glmm.fu.5yr = metaprop(rupt, total_size,
data=dat.fu.total,
studlab=paste(auth_year),
byvar = fu_5yr,
bylab = "Follow up length",
method="GLMM",
sm="PLOGIT",
method.tau = "ML",
method.ci = "WS",
pscale = 100
)
pes.summary.glmm.fu.5yr## events 95%-CI Follow up length
## Bor 2015 0.7444 [0.2535; 2.1655] Less than 5 Years
## Burns 2009 0.5780 [0.1021; 3.2011] Less than 5 Years
## Byoun 2016 1.7628 [0.9871; 3.1288] Less than 5 Years
## Choi 2018 0.5780 [0.1021; 3.2011] Greater than 5 years
## Gibbs 2004 0.0000 [0.0000; 14.8655] Greater than 5 years
## Gondar 2016 0.8152 [0.2776; 2.3691] Less than 5 Years
## Guresir 2013 0.7812 [0.2660; 2.2715] Less than 5 Years
## Irazabal 2011 0.0000 [0.0000; 7.8652] Greater than 5 years
## Jeon 2014 0.3521 [0.0966; 1.2746] Less than 5 Years
## Jiang 2013 0.0000 [0.0000; 7.1348] Less than 5 Years
## Loumiotis 2011 0.0000 [0.0000; 2.3446] Less than 5 Years
## Matsubara 2004 0.0000 [0.0000; 2.3736] Less than 5 Years
## Mizoi 1995 0.0000 [0.0000; 15.4639] Less than 5 Years
## Morita 2012 1.8658 [1.4582; 2.3845] Less than 5 Years
## Murayama 2016 2.2384 [1.6660; 3.0014] Less than 5 Years
## Oh 2013 0.0000 [0.0000; 16.8179] Less than 5 Years
## Serrone 2016 0.5155 [0.0911; 2.8615] Less than 5 Years
## So 2010 1.1450 [0.3902; 3.3118] Less than 5 Years
## Sonobe 2010 1.5625 [0.7589; 3.1897] Less than 5 Years
## Teo 2016 2.3810 [0.8130; 6.7666] Less than 5 Years
## Thien 2017 0.6173 [0.1090; 3.4133] Less than 5 Years
## Tsutsumi 2000 4.3478 [1.4896; 12.0212] Less than 5 Years
## Villablanca 2013 1.5544 [0.5300; 4.4697] Less than 5 Years
## Wiebers-R 1998 1.3503 [0.8558; 2.1244] Greater than 5 years
## Wilkinson 2018 0.0000 [0.0000; 14.8655] Greater than 5 years
## Zylkowski 2015 2.7273 [0.9318; 7.7131] Less than 5 Years
##
## Number of studies combined: k = 26
##
## events 95%-CI
## Fixed effect model 1.5091 [1.3003; 1.7509]
## Random effects model 1.1556 [0.8318; 1.6035]
##
## Quantifying heterogeneity:
## tau^2 = 0.1546; tau = 0.3932; I^2 = 45.8%; H = 1.36
##
## Quantifying residual heterogeneity:
## I^2 = 0.0% [0.0%; 40.7%]; H = 1.00 [1.00; 1.30]
##
## Test of heterogeneity:
## Q d.f. p-value Test
## 24.51 25 0.4900 Wald-type
## 46.34 25 0.0059 Likelihood-Ratio
##
## Results for subgroups (fixed effect model):
## k events 95%-CI Q I^2
## Follow up length = Less than 5 Years 21 1.5612 [1.3332; 1.8276] 22.01 51.6%
## Follow up length = Greater than 5 years 5 1.1912 [0.7611; 1.8599] 0.69 0.0%
##
## Test for subgroup differences (fixed effect model):
## Q d.f. p-value
## Between groups 1.25 1 0.2627
## Within groups 22.70 24 0.5378
##
## Results for subgroups (random effects model):
## k events 95%-CI tau^2
## Follow up length = Less than 5 Years 21 1.1638 [0.8137; 1.6619] 0.1870
## Follow up length = Greater than 5 years 5 1.1912 [0.7611; 1.8599] 0
## tau
## Follow up length = Less than 5 Years 0.4324
## Follow up length = Greater than 5 years 0
##
## Test for subgroup differences (random effects model):
## Q d.f. p-value
## Between groups 0.01 1 0.9364
##
## Details on meta-analytical method:
## - Random intercept logistic regression model
## - Maximum-likelihood estimator for tau^2
## - Logit transformation
## - Wilson Score confidence interval for individual studies
## - Continuity correction of 0.5 in studies with zero cell frequencies
## (only used to calculate individual study results)
## - Events per 100 observationsforest(pes.summary.glmm.fu.5yr,
layout = "meta",
comb.fixed = FALSE,
comb.random = TRUE,
leftlabs = c("Study", "Ruptures", "Total"),
rightcols = c("effect", "ci"),
rightlabs = c("Ruptures per 100", "95% CI"),
smlab = " ",
xlim=c(0,10),
xlab = "Rupture Proportion per 100",
pooled.events = TRUE,
) 
dat.fu.int <- dat.fu %>%
filter(fu < 60)
pes.summary.glmmint = metaprop(rupt, total_size,
data=dat.fu.int,
studlab=paste(auth_year),
method="GLMM",
sm="PLOGIT",
method.tau = "ML",
method.ci = "WS",
pscale = 100
)
pes.summary.glmmint## events 95%-CI
## Bor 2015 0.7444 [0.2535; 2.1655]
## Burns 2009 0.5780 [0.1021; 3.2011]
## Byoun 2016 1.7628 [0.9871; 3.1288]
## Gondar 2016 0.8152 [0.2776; 2.3691]
## Guresir 2013 0.7812 [0.2660; 2.2715]
## Jeon 2014 0.3521 [0.0966; 1.2746]
## Jiang 2013 0.0000 [0.0000; 7.1348]
## Loumiotis 2011 0.0000 [0.0000; 2.3446]
## Matsubara 2004 0.0000 [0.0000; 2.3736]
## Mizoi 1995 0.0000 [0.0000; 15.4639]
## Morita 2012 1.8658 [1.4582; 2.3845]
## Murayama 2016 2.2384 [1.6660; 3.0014]
## Oh 2013 0.0000 [0.0000; 16.8179]
## Serrone 2016 0.5155 [0.0911; 2.8615]
## So 2010 1.1450 [0.3902; 3.3118]
## Sonobe 2010 1.5625 [0.7589; 3.1897]
## Teo 2016 2.3810 [0.8130; 6.7666]
## Thien 2017 0.6173 [0.1090; 3.4133]
## Tsutsumi 2000 4.3478 [1.4896; 12.0212]
## Villablanca 2013 1.5544 [0.5300; 4.4697]
## Zylkowski 2015 2.7273 [0.9318; 7.7131]
##
## Number of studies combined: k = 21
##
## events 95%-CI
## Fixed effect model 1.5612 [1.3332; 1.8276]
## Random effects model 1.1638 [0.8137; 1.6619]
##
## Quantifying heterogeneity:
## tau^2 = 0.1870; tau = 0.4324; I^2 = 51.6%; H = 1.44
##
## Test of heterogeneity:
## Q d.f. p-value Test
## 22.01 20 0.3402 Wald-type
## 41.91 20 0.0028 Likelihood-Ratio
##
## Details on meta-analytical method:
## - Random intercept logistic regression model
## - Maximum-likelihood estimator for tau^2
## - Logit transformation
## - Wilson Score confidence interval for individual studies
## - Continuity correction of 0.5 in studies with zero cell frequencies
## (only used to calculate individual study results)
## - Events per 100 observationsSubgroup forest plots for patients with less than 5 years of follow up
forest(pes.summary.glmmint,
layout = "meta",
comb.fixed = FALSE,
comb.random = TRUE,
leftlabs = c("Study", "Ruptures", "Total"),
rightcols = c("effect", "ci"),
rightlabs = c("Ruptures per 100", "95% CI"),
smlab = " ",
xlim=c(0,10),
xlab = "Rupture Proportion per 100",
pooled.events = TRUE,
) 
Meta-analysis of proportions for patients with 5 years or more of follow up
dat.fu.long <- dat.fu %>%
filter(fu > 60)
pes.summary.glmmlong = metaprop(rupt, total_size,
data=dat.fu.long,
studlab=paste(auth_year),
method="GLMM",
sm="PLOGIT",
method.tau = "ML",
method.ci = "WS",
pscale = 100
)
pes.summary.glmmlong## events 95%-CI
## Choi 2018 0.5780 [0.1021; 3.2011]
## Gibbs 2004 0.0000 [0.0000; 14.8655]
## Irazabal 2011 0.0000 [0.0000; 7.8652]
## Wiebers-R 1998 1.3503 [0.8558; 2.1244]
## Wilkinson 2018 0.0000 [0.0000; 14.8655]
##
## Number of studies combined: k = 5
##
## events 95%-CI
## Fixed effect model 1.1912 [0.7611; 1.8599]
## Random effects model 1.1912 [0.7611; 1.8599]
##
## Quantifying heterogeneity:
## tau^2 = 0; tau = 0; I^2 = 0.0%; H = 1.00
##
## Test of heterogeneity:
## Q d.f. p-value Test
## 0.69 4 0.9525 Wald-type
## 3.09 4 0.5429 Likelihood-Ratio
##
## Details on meta-analytical method:
## - Random intercept logistic regression model
## - Maximum-likelihood estimator for tau^2
## - Logit transformation
## - Wilson Score confidence interval for individual studies
## - Continuity correction of 0.5 in studies with zero cell frequencies
## (only used to calculate individual study results)
## - Events per 100 observationsforest(pes.summary.glmmlong,
layout = "meta",
comb.fixed = FALSE,
comb.random = TRUE,
leftlabs = c("Study", "Ruptures", "Total"),
rightcols = c("effect", "ci"),
rightlabs = c("Ruptures per 100", "95% CI"),
smlab = " ",
xlim=c(0,10),
xlab = "Rupture Proportion per 100",
pooled.events = TRUE,
) 
The data for risk of rupture beyond 5 years is extremely sparse, and quite imprecise. Pooling is not considered appropriate, and forest plots are not required. The results will be described qualitatively.
It is noted however, that the result in the long term follow up subgroup may be due to mainly patients with prior history of SAH and long term follow up. Explore this by examining only patients with prior SAH and long term follow up.
dat.psah.fu <- sizedata %>%
mutate(fu = coalesce(fu_mean_tot,fu_med_tot)) %>%
unite(auth_year, c(auth, pub), sep = " ", remove = FALSE) %>%
mutate(total_anr = coalesce(num_anr,num)) %>%
drop_na(total_anr) %>%
rename(ai = rupt_psah) %>%
rename(rupt.total = rupt) %>%
mutate(ci = rupt.total - ai) %>%
rename(n1i.temp = psah) %>%
mutate(prop = psah_tot / num_tot) %>%
mutate(num_anr_psah = prop * total_anr) %>%
mutate(n1i = coalesce(n1i.temp,num_anr_psah)) %>%
mutate(n2i = total_anr - n1i) %>%
mutate(bi = n1i - ai) %>%
mutate(di = n2i - ci) %>%
select(auth_year, ai, bi, ci, di, n1i, n2i, fu) %>%
drop_na(ai, bi, ci, di) %>%
mutate_if(is.numeric, round, 0)How about patients with intermediate follow up?
dat.psah.fu.int <- dat.psah.fu %>%
filter(fu < 60)How about patients with long term follow up?
dat.psah.fu.long <- dat.psah.fu %>%
filter(fu > 60)Examine dat.psah.fu.long and review ai and ci.
ci = patients without prior history of SAH who were included in studies with greater than 5 years of follow up who susbsequently ruptured.
See that only 2 patients without history of prior SAH ruptured on long term follow up. Although the data are sparse, there is minimal extractable data on the outcomes of unruptured aneurysms beyond 5 years.
Meta-regression
Meta-regression
- Prespecified metagression to consider whether proportion of patients with exposure to prior SAH in the studies was an important source of heterogeneity.
- Prespecified metagression to consider whether proportion of patients with aneurysms 5mm and less in the studies was an important source of heterogeneity.
- Prespecified metagression to consider whether proportion of patients with aneurysms 3mm and less in the studies was an important source of heterogeneity.
- Post hoc metagression to consider whether duration of follow up in the studies analysed as a continuous variable was an important source of heterogeneity.
- Post hoc metagression to consider whether year of last enrollment in the studies analysed was an important source of heterogeneity.
To reduce the risk of false positives, and to ensure that there are adequate studies to consider pos hoc meta-regression, only one posthoc potential effect modifier is considered for metaregression each time.
In the multivariable meta-regression model only 2 potential effect modifiers are considered at a time.
Regardless, the outcomes of these subgroup analysis or meta-regression analyses are limited. They are entirely observational, susceptible to confounding bias from other study level characteristics, the level of statistical significance of the results within subgroups cannot be compared, and may not be generalisable since there is aggregation bias due to study-level data being used for analysis and not patient level data.
Exploratory meta-regression to explore length of follow up on rupture outcomes, including all studies
dat.fu.metareg <- dat %>%
mutate(fu = coalesce(fu_mean_tot,fu_med_tot)/12) %>%
drop_na(total_size) %>%
unite(auth_year, c(auth, pub), sep = " ", remove = FALSE) %>%
select(auth_year, rupt, total_size, fu) %>%
drop_na(fu)fu.metareg = metaprop(rupt, total_size,
data=dat.fu.metareg,
studlab=paste(auth_year),
method="GLMM",
sm="PLOGIT",
method.tau = "ML",
method.ci = "WS",
pscale = 100
)fu.metareg.result <- metareg(fu.metareg, ~fu)
fu.metareg.result##
## Mixed-Effects Model (k = 26; tau^2 estimator: ML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.1553
## tau (square root of estimated tau^2 value): 0.3940
## I^2 (residual heterogeneity / unaccounted variability): 41.7685%
## H^2 (unaccounted variability / sampling variability): 1.7173
##
## Tests for Residual Heterogeneity:
## Wld(df = 24) = 23.8497, p-val = 0.4702
## LRT(df = 24) = 45.8505, p-val = 0.0046
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 0.0031, p-val = 0.9554
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt -4.4631 0.3066 -14.5574 <.0001 -5.0640 -3.8622 ***
## fu 0.0040 0.0717 0.0559 0.9554 -0.1364 0.1445
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1bubble(fu.metareg.result, studlab=dat.fu.metareg$auth_year, xlim=c(0, 10), main="Meta-regression plot for length of FU")
Overall, the metaregression result reveals that there is no association between study follow up length and rupture risk.
estimate 0.0040 (95% CI -0.1364 to 0.1445) p value <.0001
Metaregression to explore exposure to prior SAH
We can analyse whether there is a relationship between the proportion of patients with prior SAH included in the study (potential effect modifier), and the the outcome (rupture proportion).
The choice of this potential effect modifier is:
- Pre-specified
- Biologically plausible
- Supported by published observational data
- Relevant to physicians
Since we are analysing at a study level, and we do not know the characteristics of individual patients, this is an exploratory and hypothesis generating analysis.
Note that
ai = prior SAH and rupture +ve = variable rupt_psah bi = prior SAH and rupture -ve ci = no prior SAH and rupture +ve = rupt minus rupt_psah di = no prior SAH, and rupture -ve
n1i = ai + bi = total aneurysms with prior SAH n2i = ci + di = total aneurysms without prior SAH
n1i + n2i = total_anr included with and without prior SAH
total number of events = total number of ruptures = variable rupt
dat.psah.metareg <- dat %>%
unite(auth_year, c(auth, pub), sep = " ", remove = FALSE) %>%
select(auth_year, rupt, total_size, prop_psah) %>%
drop_na(total_size) %>%
drop_na(prop_psah)psah.metareg = metaprop(rupt, total_size,
data=dat.psah.metareg,
studlab=paste(auth_year),
method="GLMM",
sm="PLOGIT",
method.tau = "ML",
method.ci = "WS",
pscale = 100
)
psah.metareg.result <- metareg(psah.metareg, ~prop_psah)
psah.metareg.result##
## Mixed-Effects Model (k = 27; tau^2 estimator: ML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.2152
## tau (square root of estimated tau^2 value): 0.4639
## I^2 (residual heterogeneity / unaccounted variability): 49.0663%
## H^2 (unaccounted variability / sampling variability): 1.9633
##
## Tests for Residual Heterogeneity:
## Wld(df = 25) = 26.1566, p-val = 0.3993
## LRT(df = 25) = 49.4678, p-val = 0.0025
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 0.8631, p-val = 0.3529
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt -4.5639 0.2395 -19.0593 <.0001 -5.0332 -4.0946 ***
## prop_psah 0.9979 1.0742 0.9290 0.3529 -1.1074 3.1032
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1bubble(psah.metareg.result, studlab=dat.psah.metareg$auth_year, xlim=c(0, 1), main="Meta-regression plot for prior SAH")
Here again, metaregression reveals that the proportion of patients with prior SAH included in the study (potential effect modifier) is not associated with the outcome (rupture proportion).
Metaregression to explore size categorised at 5mm
We can analyse whether there is a relationship between the proportion of patients with aneurysms measuring 5mm and less included in the study (potential effect modifier), and the the outcome (rupture proportion).
The choice of this potential effect modifier is:
- Pre-specified
- Biologically plausible
- Supported by published observational data
- Relevant to physicians
Since we are analysing at a study level, and we do not know the aneurysm characteristics of individual patients, this is an exploratory and hypothesis generating analysis.
sizedata10 <- filter(maindata, size == 10)
dat10.metareg <- sizedata10 %>%
slice(-12) %>%
mutate(total_size10 = coalesce(num_anr,num)) %>%
drop_na(total_size10) %>%
unite(auth_year, c(auth, pub), sep = " ", remove = FALSE) %>%
select(auth_year, rupt, total_size10) %>%
rename(rupt10 = rupt)
dat5.metareg <- dat5 %>%
rename(rupt5 = rupt) %>%
rename(total_size5 = total_size)
dat.size5.metareg <- left_join(dat5.metareg, dat10.metareg) %>%
mutate(prop_size5 = total_size5 / total_size10)## Joining, by = "auth_year"size5.metareg = metaprop(rupt10, total_size10,
data=dat.size5.metareg,
studlab=paste(auth_year),
method="GLMM",
sm="PLOGIT",
method.tau = "ML",
method.ci = "WS",
pscale = 100
)
size5mm.metareg.result <- metareg(size5.metareg, ~prop_size5)
size5mm.metareg.result##
## Mixed-Effects Model (k = 22; tau^2 estimator: ML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.2385
## tau (square root of estimated tau^2 value): 0.4883
## I^2 (residual heterogeneity / unaccounted variability): 45.9817%
## H^2 (unaccounted variability / sampling variability): 1.8512
##
## Tests for Residual Heterogeneity:
## Wld(df = 20) = 23.6374, p-val = 0.2586
## LRT(df = 20) = 41.6291, p-val = 0.0031
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 0.3404, p-val = 0.5596
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt -4.8648 0.8567 -5.6787 <.0001 -6.5439 -3.1858 ***
## prop_size5 0.6023 1.0324 0.5834 0.5596 -1.4211 2.6258
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1bubble(size5mm.metareg.result, studlab=dat.size5.metareg$auth_year, xlim=c(0, 1), main="Meta-regression plot for 5mm Aneurysm size")
There is no evidence of that there is an association between pooled rupture risk and theproportion of patients with aneursysms 5mm and less included in the studies
Metaregression to explore size categorised at 3mm
sizedata10 <- filter(maindata, size == 10)
dat10.metareg <- sizedata10 %>%
slice(-12) %>%
mutate(total_size10 = coalesce(num_anr,num)) %>%
drop_na(total_size10) %>%
unite(auth_year, c(auth, pub), sep = " ", remove = FALSE) %>%
select(auth_year, rupt, total_size10) %>%
rename(rupt10 = rupt)
dat3.metareg <- dat3 %>%
rename(rupt3 = rupt) %>%
rename(total_size3 = total_size)
dat.size3.metareg <- left_join(dat3.metareg, dat10.metareg) %>%
mutate(prop_size3 = total_size3 / total_size10)## Joining, by = "auth_year"size3.metareg = metaprop(rupt10, total_size10,
data=dat.size3.metareg,
studlab=paste(auth_year),
method="GLMM",
sm="PLOGIT",
method.tau = "ML",
method.ci = "WS",
pscale = 100
)
size3mm.metareg.result <- metareg(size3.metareg, ~prop_size3)
size3mm.metareg.result##
## Mixed-Effects Model (k = 15; tau^2 estimator: ML)
##
## tau^2 (estimated amount of residual heterogeneity): 0
## tau (square root of estimated tau^2 value): 0
## I^2 (residual heterogeneity / unaccounted variability): 0.0000%
## H^2 (unaccounted variability / sampling variability): 1.0000
##
## Tests for Residual Heterogeneity:
## Wld(df = 13) = 5.7603, p-val = 0.9544
## LRT(df = 13) = 10.5569, p-val = 0.6479
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 0.3458, p-val = 0.5565
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt -4.7066 0.3961 -11.8822 <.0001 -5.4829 -3.9302 ***
## prop_size3 0.5817 0.9893 0.5880 0.5565 -1.3572 2.5207
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1bubble(size3mm.metareg.result, studlab=dat.size3.metareg$auth_year, xlim=c(0, 1), main="Meta-regression plot for 3mm Aneurysm size")
There is no evidence that the pooled rupture risk varied by the proportion of patients with aneursysms 3mm and less included in the studies.
Metaregression to explore size end of study enrollment
We can analyse whether there is a relationship between the year the last patient was included in the study (potential effect modifier), and the the outcome (rupture proportion).
The choice of this potential effect modifier is:
- Exploratory
- Biologically plausible - since comorbidity profiles, and medical management of comorbidities has changed
- Relevant to physicians
Since we are analysing at a study level, this is an exploratory and hypothesis generating analysis.
dat.year.metareg <- sizedata %>%
slice(-12) %>%
mutate(total_size = coalesce(num_anr,num)) %>%
drop_na(total_size) %>%
unite(auth_year, c(auth, pub), sep = " ", remove = FALSE) %>%
select(auth_year, end, rupt, total_size) %>%
drop_na(end)year.metareg = metaprop(rupt, total_size,
data=dat.year.metareg,
studlab=paste(auth_year),
method="GLMM",
sm="PLOGIT",
method.tau = "ML",
method.ci = "WS",
pscale = 100
)
year.metareg.result <- metareg(year.metareg, ~end)
year.metareg.result##
## Mixed-Effects Model (k = 25; tau^2 estimator: ML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.1377
## tau (square root of estimated tau^2 value): 0.3711
## I^2 (residual heterogeneity / unaccounted variability): 41.0342%
## H^2 (unaccounted variability / sampling variability): 1.6959
##
## Tests for Residual Heterogeneity:
## Wld(df = 23) = 24.4260, p-val = 0.3805
## LRT(df = 23) = 40.8728, p-val = 0.0122
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 2.6952, p-val = 0.1007
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 71.4488 46.1113 1.5495 0.1213 -18.9278 161.8253
## end -0.0378 0.0230 -1.6417 0.1007 -0.0828 0.0073
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1```{r metaregression - option 2 - bubble plot - year} 2
bubble(year.metareg.result, studlab=dat.year.metareg$auth_year, xlim=c(1980, 2020), main=“Meta-regression plot for Study End Year”)
This plot suggests that there is no association between end of study year and rupture risk of aneurysms.
We can repeat this excluding Juvela.
```r
dat.year.metareg.juvela <- sizedata %>%
slice(-12) %>%
mutate(total_size = coalesce(num_anr,num)) %>%
drop_na(total_size) %>%
unite(auth_year, c(auth, pub), sep = " ", remove = FALSE) %>%
select(auth_year, end, rupt, total_size) %>%
drop_na(end)year.metareg.juvela = metaprop(rupt, total_size,
data=dat.year.metareg.juvela,
studlab=paste(auth_year),
method="GLMM",
sm="PLOGIT",
method.tau = "ML",
method.ci = "WS",
pscale = 100
)
year.metareg.result.juvela <- metareg(year.metareg.juvela, ~end)
year.metareg.result.juvela##
## Mixed-Effects Model (k = 25; tau^2 estimator: ML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.1377
## tau (square root of estimated tau^2 value): 0.3711
## I^2 (residual heterogeneity / unaccounted variability): 41.0342%
## H^2 (unaccounted variability / sampling variability): 1.6959
##
## Tests for Residual Heterogeneity:
## Wld(df = 23) = 24.4260, p-val = 0.3805
## LRT(df = 23) = 40.8728, p-val = 0.0122
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 2.6952, p-val = 0.1007
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 71.4488 46.1113 1.5495 0.1213 -18.9278 161.8253
## end -0.0378 0.0230 -1.6417 0.1007 -0.0828 0.0073
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1bubble(year.metareg.result.juvela, studlab=dat.year.metareg.juvela$auth_year, xlim=c(1985, 2015), main="Meta-regression plot for Study End Year excluding Juvela")
Once Juvela is excluded, then there is no evidence of an association between rupture risk and study end year.
Metaregression to explore Japanese population vs non-Japanese
We can analyse whether there is a relationship between whether the source population for the study was Japan (potential effect modifier), and the outcome (rupture proportion).
The choice of this potential effect modifier is:
- Exploratory
- Biologically plausible - since it remains unclear why Japanese population has a higher SAH rate compared to the rest of the world
- Relevant to physicians
Since we are analysing at a study level, this is an exploratory and hypothesis generating analysis.
dat.jap.metareg <- dat %>%
drop_na(total_size) %>%
unite(auth_year, c(auth, pub), sep = " ", remove = FALSE) %>%
select(auth_year, pop, rupt, total_size) %>%
drop_na(pop)jap.metareg = metaprop(rupt, total_size,
data=dat.jap.metareg,
studlab=paste(auth_year),
method="GLMM",
sm="PLOGIT",
method.tau = "ML",
method.ci = "WS",
pscale = 100
)
jap.metareg.result <- metareg(jap.metareg, ~pop)
jap.metareg.result##
## Mixed-Effects Model (k = 29; tau^2 estimator: ML)
##
## tau^2 (estimated amount of residual heterogeneity): 0
## tau (square root of estimated tau^2 value): 0
## I^2 (residual heterogeneity / unaccounted variability): 0.0000%
## H^2 (unaccounted variability / sampling variability): 1.0000
##
## Tests for Residual Heterogeneity:
## Wld(df = 27) = 20.8937, p-val = 0.7912
## LRT(df = 27) = 36.3088, p-val = 0.1087
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 12.9665, p-val = 0.0003
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt -4.4838 0.1267 -35.3899 <.0001 -4.7321 -4.2355 ***
## popJapanese 0.5659 0.1571 3.6009 0.0003 0.2579 0.8739 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1bubble(jap.metareg.result, studlab=dat.jap.metareg$auth_year, xlim=c(1960, 2020), main="Meta-regression plot for Study End Year")
Multivariable meta-regression
In the previous scenarios we have considered the impact of a single study-level effect modifier on the outcome of rupture risk.
Given that we have identified that when all studies are included, the proportion of patients with prior SAH exposure and year of study end are associated with rupture risk, we can examine this in a multivariable meta-regression model.
dat.psah.year.metareg <- sizedata %>%
unite(auth_year, c(auth, pub), sep = " ", remove = FALSE) %>%
mutate(total_anr = coalesce(num_anr,num)) %>%
drop_na(total_anr) %>%
rename(ai = rupt_psah) %>%
rename(rupt.total = rupt) %>%
mutate(ci = rupt.total - ai) %>%
rename(n1i.temp = psah) %>%
mutate(prop = psah_tot / num_tot) %>%
mutate(num_anr_psah = prop * total_anr) %>%
mutate(n1i = coalesce(n1i.temp,num_anr_psah)) %>%
mutate(n2i = total_anr - n1i) %>%
mutate(bi = n1i - ai) %>%
mutate(di = n2i - ci) %>%
select(auth_year, rupt.total, total_anr, prop, end) psah.year.metareg = metaprop(rupt.total, total_anr,
data=dat.psah.year.metareg,
studlab=paste(auth_year),
method="GLMM",
sm="PLOGIT",
method.tau = "ML",
method.ci = "WS",
pscale = 100
)
psah.year.metareg.result <- metareg(psah.year.metareg, ~ prop + end)## Warning in rma.glmm(xi = event[!exclude], ni = n[!exclude], data = dataset, :
## Studies with NAs omitted from model fitting.## Warning in rma.glmm(xi = event[!exclude], ni = n[!exclude], data = dataset, :
## Some yi/vi values are NA.psah.year.metareg.result##
## Mixed-Effects Model (k = 24; tau^2 estimator: ML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.3194
## tau (square root of estimated tau^2 value): 0.5651
## I^2 (residual heterogeneity / unaccounted variability): 58.7915%
## H^2 (unaccounted variability / sampling variability): 2.4267
##
## Tests for Residual Heterogeneity:
## Wld(df = 21) = 58.5954, p-val < .0001
## LRT(df = 21) = 79.7854, p-val < .0001
##
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 25.6939, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 106.0048 70.7508 1.4983 0.1341 -32.6642 244.6739
## prop 1.4246 1.3481 1.0567 0.2906 -1.2176 4.0668
## end -0.0551 0.0352 -1.5654 0.1175 -0.1241 0.0139
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1bubble(psah.year.metareg.result, studlab=dat.psah.year.metareg$auth_year, xlim=c(0, 1), main="Multivariable Meta-regression plot for prior SAH and Study End")## Warning in bubble.metareg(psah.year.metareg.result, studlab =
## dat.psah.year.metareg$auth_year, : Only first covariate in meta-regression
## ('prop') considered in bubble plot. No regression line plotted.
Proportion of prior SAH patients in each study and year of study end are not associated with rupture risk when considered in a multivariable meta-regression analysis.
consider additional multivariable meta-regression with proportion of prior SAH and study follow up
Sensitivity analysis
Sensitivity analysis is used to assess the robustness of the result, and helps strengthed or weaken the conclusions that can be drawn from the meta-analysis.
If sensitivity analyses show that the overall result and conclusions are not affected by the different decisions that could be made during the review process, the results of the review can be regarded with a higher degree of certainty.
Similarly, if sensitivity analyses show that the overall result and conclusions are affected, the results should be interpreted with caution, and can be used to generate hypothesis for additional research.
Sensitivity analysis differs from subgroup analysis in that there is no assessment of the treatment effect in the removed studies, nor is there formal statistical assessment across the included and excluded studies. Instead, there is an informal comparison by recalculating the pooled effect size.
Sensitivity analysis
- Post hoc sensitivity analysis to consider whether the pooled rupture risk varies with utilisation of a fixed or random effects meta-analysis method.
- Post hoc sensitivity analysis to consider whether the pooled rupture risk varies with exclusion of studies rated as poor or fair quality using AHRQ thresholds of NOS quality scoring.
- Post hoc sensitivity analysis to consider whether the pooled rupture risk varies with exclusion of studies with sample sizes of less than 500 aneurysms.
- Post hoc sensitivity analysis to consider whether the pooled rupture risk varies with exclusion of studie that are identified as contributing considerably to the overall heterogeneity and having a strong influence on the overall result using a Baujat plot.
- Repeat metagression with exlusion of Juvela et al to consider whether proportion of patients with exposure to prior SAH in the study was an important source of heterogeneity.
- Repeat metagression with exlusion of Juvela et al to consider whether duration of follow up in the study analysed as a continuous variable was an important source of heterogeneity
- Repeat metagression with exlusion of Juvela et al to consider whether year of last enrollment in the study analysed as a continuous variable was an important source of heterogeneity
- Post hoc multivariable metaregression to consider whether inclusion of 2 covariates - proportion of patients with exposure to prior SAH in the study and year of last enrollment in the study - was associated with pooled rupture risk and a source of heterogeniety.
Sensitivity analysis - consider fixed vs random effects meta-analysis decision
If there is heterogeniety, the random effect meta-analysis weights studies more equally. However, if there are marked small study effects, and only 1 high quality study, then this may not be reasonable. To synthesise the data, the author must make a decision whether this is appropriate. This can be done by comparing the results of the random and fixed effects meta-analysis.
pes.summary.glmm## events 95%-CI
## Bor 2015 0.7444 [0.2535; 2.1655]
## Broderick 2009 1.8519 [0.5093; 6.5019]
## Burns 2009 0.5780 [0.1021; 3.2011]
## Byoun 2016 1.7628 [0.9871; 3.1288]
## Choi 2018 0.5780 [0.1021; 3.2011]
## Gibbs 2004 0.0000 [0.0000; 14.8655]
## Gondar 2016 0.8152 [0.2776; 2.3691]
## Guresir 2013 0.7812 [0.2660; 2.2715]
## Irazabal 2011 0.0000 [0.0000; 7.8652]
## Jeon 2014 0.3521 [0.0966; 1.2746]
## Jiang 2013 0.0000 [0.0000; 7.1348]
## Loumiotis 2011 0.0000 [0.0000; 2.3446]
## Matsubara 2004 0.0000 [0.0000; 2.3736]
## Matsumoto 2013 2.6549 [0.9069; 7.5160]
## Mizoi 1995 0.0000 [0.0000; 15.4639]
## Morita 2012 1.8658 [1.4582; 2.3845]
## Murayama 2016 2.2384 [1.6660; 3.0014]
## Oh 2013 0.0000 [0.0000; 16.8179]
## Serrone 2016 0.5155 [0.0911; 2.8615]
## So 2010 1.1450 [0.3902; 3.3118]
## Sonobe 2010 1.5625 [0.7589; 3.1897]
## Teo 2016 2.3810 [0.8130; 6.7666]
## Thien 2017 0.6173 [0.1090; 3.4133]
## Tsukahara 2005 3.4722 [1.4921; 7.8703]
## Tsutsumi 2000 4.3478 [1.4896; 12.0212]
## Villablanca 2013 1.5544 [0.5300; 4.4697]
## Wiebers-R 1998 1.3503 [0.8558; 2.1244]
## Wilkinson 2018 0.0000 [0.0000; 14.8655]
## Zylkowski 2015 2.7273 [0.9318; 7.7131]
##
## Number of studies combined: k = 29
##
## events 95%-CI
## Fixed effect model 1.5475 [1.3390; 1.7879]
## Random effects model 1.2625 [0.9445; 1.6858]
##
## Quantifying heterogeneity:
## tau^2 = 0.1359; tau = 0.3687; I^2 = 41.8%; H = 1.31
##
## Test of heterogeneity:
## Q d.f. p-value Test
## 27.45 28 0.4940 Wald-type
## 49.86 28 0.0067 Likelihood-Ratio
##
## Details on meta-analytical method:
## - Random intercept logistic regression model
## - Maximum-likelihood estimator for tau^2
## - Logit transformation
## - Wilson Score confidence interval for individual studies
## - Continuity correction of 0.5 in studies with zero cell frequencies
## (only used to calculate individual study results)
## - Events per 100 observationsHere the results are not clinically significant with the overall same rupture proportion of 1-2% with either method. Moreover, in our case, a random effects meta-analysis is more appropriate due to existence of clinical and methodological heterogeneity. Regardless, other sources of heterogeneity still need to be explored.
Sensitivity analysis - exclude based on study quality
Quality assessment was performed using the Newastle Ottowa Scale as recommended by Cochrane; this is also the most commonly utilised for observational studies.
This can be converted to the The Agency for Healthcare Research and Quality within the United States Department of Health and Human Services (AHRQ) standards using the following thresholds.
Good quality: 3 or 4 stars in selection domain AND 1 or 2 stars in comparability domain AND 2 or 3 stars in outcome/exposure domain
Fair quality: 2 stars in selection domain AND 1 or 2 stars in comparability domain AND 2 or 3 stars in outcome/exposure domain
Poor quality: 0 or 1 star in selection domain OR 0 stars in comparability domain OR 0 or 1 stars in outcome/exposure domain
dat.sens.nos <- sizedata %>%
filter(nos_select == 3 | nos_select == 4 ) %>%
filter(nos_compare == 1 | nos_compare == 2) %>%
filter(nos_outcome == 2 | nos_outcome == 3) %>%
mutate(ahrq = "Good")
dat.sens.nos <- dat.sens.nos %>%
mutate(total_size = coalesce(num_anr,num)) %>%
drop_na(total_size) %>%
unite(auth_year, c(auth, pub), sep = " ", remove = FALSE) %>%
select(auth_year, rupt, total_size)pes.sens.nos = metaprop(rupt, total_size,
data=dat.sens.nos,
studlab=paste(auth_year),
method="GLMM",
sm="PLOGIT",
method.tau = "ML",
method.ci = "WS",
pscale = 100
)
pes.sens.nos## events 95%-CI
## Irazabal 2011 0.0000 [ 0.0000; 7.8652]
## Juvela 2013 22.9630 [16.6751; 30.7469]
## Matsumoto 2013 2.6549 [ 0.9069; 7.5160]
## Mizoi 1995 0.0000 [ 0.0000; 15.4639]
## Morita 2012 1.8658 [ 1.4582; 2.3845]
## Murayama 2016 2.2384 [ 1.6660; 3.0014]
## Serrone 2016 0.5155 [ 0.0911; 2.8615]
## So 2010 1.1450 [ 0.3902; 3.3118]
## Sonobe 2010 1.5625 [ 0.7589; 3.1897]
## Wiebers-R 1998 1.6901 [ 1.0717; 2.6558]
##
## Number of studies combined: k = 10
##
## events 95%-CI
## Fixed effect model 2.2320 [1.9215; 2.5912]
## Random effects model 1.9027 [0.8747; 4.0888]
##
## Quantifying heterogeneity:
## tau^2 = 1.1526; tau = 1.0736; I^2 = 93.8%; H = 4.03
##
## Test of heterogeneity:
## Q d.f. p-value Test
## 155.62 9 < 0.0001 Wald-type
## 108.28 9 < 0.0001 Likelihood-Ratio
##
## Details on meta-analytical method:
## - Random intercept logistic regression model
## - Maximum-likelihood estimator for tau^2
## - Logit transformation
## - Wilson Score confidence interval for individual studies
## - Continuity correction of 0.5 in studies with zero cell frequencies
## (only used to calculate individual study results)
## - Events per 100 observationsHere we can see the point estimate and the confidence intervals varies from
All studies 1.2198 [0.7739; 1.9177] Good Quality studies 1.9027 [0.8747; 4.0888]
There is even higher heterogeniety, and thus results of the original pooled analysis need to be interpreted with caution.
Sensitivity analysis - exclude based on small study samples
dat.sens.samplesize <- sizedata %>%
mutate(total_size = coalesce(num_anr,num)) %>%
drop_na(total_size) %>%
unite(auth_year, c(auth, pub), sep = " ", remove = FALSE) %>%
filter(total_size > 500) %>%
select(auth_year, rupt, total_size)pes.sens.samplesize = metaprop(rupt, total_size,
data=dat.sens.samplesize,
studlab=paste(auth_year),
method="GLMM",
sm="PLOGIT",
method.tau = "ML",
method.ci = "WS",
pscale = 100
)
pes.sens.samplesize## events 95%-CI
## Byoun 2016 1.7628 [0.9871; 3.1288]
## Jeon 2014 0.3521 [0.0966; 1.2746]
## Morita 2012 1.8658 [1.4582; 2.3845]
## Murayama 2016 2.2384 [1.6660; 3.0014]
## Wiebers-R 1998 1.6901 [1.0717; 2.6558]
##
## Number of studies combined: k = 5
##
## events 95%-CI
## Fixed effect model 1.8131 [1.5346; 2.1411]
## Random effects model 1.7312 [1.2165; 2.4584]
##
## Quantifying heterogeneity:
## tau^2 = 0.0250; tau = 0.1582; I^2 = 35.8%; H = 1.25
##
## Test of heterogeneity:
## Q d.f. p-value Test
## 7.28 4 0.1218 Wald-type
## 12.13 4 0.0164 Likelihood-Ratio
##
## Details on meta-analytical method:
## - Random intercept logistic regression model
## - Maximum-likelihood estimator for tau^2
## - Logit transformation
## - Wilson Score confidence interval for individual studies
## - Events per 100 observationsHere we can see the point estimate and the confidence intervals varies from
All studies 1.2198 [0.7739; 1.9177] Studies with sample size > 500 aneurysms 1.7312 [1.2165; 2.4584]
Notably, the confidence intervals overlap, and I2 has reduced. Overall, suggests that the original pooled result needs to be interpreted with caution.
Sensitivity analysis - using Baujat plots to identify which studies contribute most to heterogeniety
These plots shows the contribution of each study to the overall Q-test statistic for heterogeneity on the horizontal axis versus the influence of each study (defined as the standardized squared difference between the overall estimate based on a fixed-effects model with and without the study included in the model) on the vertical axis. Baujat 2002
baujat(pes.summary.glmm)
This shows that Juvela et al is the major source of between study heterogeneity.
We can use these sources of heterogeniety to assess for moderating variables that may contribute to the heterogeneity.
Re-run analysis excluding Juvela - post hoc sensitivity analysis
dat.juvela <- slice(dat, -12)
pes.summary.glmm.juvela = metaprop(rupt, total_size,
data=dat.juvela,
studlab=paste(auth_year),
sm="PLOGIT",
method.tau = "ML",
method.ci = "WS",
pscale = 100
)
pes.summary.glmm.juvela## events 95%-CI
## Bor 2015 0.7444 [0.2535; 2.1655]
## Broderick 2009 1.8519 [0.5093; 6.5019]
## Burns 2009 0.5780 [0.1021; 3.2011]
## Byoun 2016 1.7628 [0.9871; 3.1288]
## Choi 2018 0.5780 [0.1021; 3.2011]
## Gibbs 2004 0.0000 [0.0000; 14.8655]
## Gondar 2016 0.8152 [0.2776; 2.3691]
## Guresir 2013 0.7812 [0.2660; 2.2715]
## Irazabal 2011 0.0000 [0.0000; 7.8652]
## Jeon 2014 0.3521 [0.0966; 1.2746]
## Jiang 2013 0.0000 [0.0000; 7.1348]
## Matsubara 2004 0.0000 [0.0000; 2.3736]
## Matsumoto 2013 2.6549 [0.9069; 7.5160]
## Mizoi 1995 0.0000 [0.0000; 15.4639]
## Morita 2012 1.8658 [1.4582; 2.3845]
## Murayama 2016 2.2384 [1.6660; 3.0014]
## Oh 2013 0.0000 [0.0000; 16.8179]
## Serrone 2016 0.5155 [0.0911; 2.8615]
## So 2010 1.1450 [0.3902; 3.3118]
## Sonobe 2010 1.5625 [0.7589; 3.1897]
## Teo 2016 2.3810 [0.8130; 6.7666]
## Thien 2017 0.6173 [0.1090; 3.4133]
## Tsukahara 2005 3.4722 [1.4921; 7.8703]
## Tsutsumi 2000 4.3478 [1.4896; 12.0212]
## Villablanca 2013 1.5544 [0.5300; 4.4697]
## Wiebers-R 1998 1.3503 [0.8558; 2.1244]
## Wilkinson 2018 0.0000 [0.0000; 14.8655]
## Zylkowski 2015 2.7273 [0.9318; 7.7131]
##
## Number of studies combined: k = 28
##
## events 95%-CI
## Fixed effect model 1.5690 [1.3576; 1.8127]
## Random effects model 1.3184 [1.0002; 1.7359]
##
## Quantifying heterogeneity:
## tau^2 = 0.1125; tau = 0.3354; I^2 = 38.1%; H = 1.27
##
## Test of heterogeneity:
## Q d.f. p-value Test
## 27.45 27 0.4399 Wald-type
## 44.83 27 0.0169 Likelihood-Ratio
##
## Details on meta-analytical method:
## - Random intercept logistic regression model
## - Maximum-likelihood estimator for tau^2
## - Logit transformation
## - Wilson Score confidence interval for individual studies
## - Continuity correction of 0.5 in studies with zero cell frequencies
## (only used to calculate individual study results)
## - Events per 100 observationsThis exploratory analaysis demonstrates greater homogeniety, with reduced and now moderate I2 and higher Q-statistic.
The point estimate and confidence intervals have changed slightly, which confirms the influence of the Juvela study.
Including Juvela - 1.2198 (0.7739; 1.9177) Excluding Juvela - 1.2849 (0.9547; 1.7274)
Nonetheless, the change is not significant from a clinical application point of view, with overall rupture risk still 1-2% overall.
Rerun Baujat plots excluing Juvela - exploratory
baujat(pes.summary.glmm.juvela,
xlim=c(0,40),
ylim=c(0,40)
)
Keeping the influence axis the same, the significant improvement in heterogeniety is noted. This can be re-processed with smaller scales to more closely explore the result.
baujat(pes.summary.glmm.juvela,
xlim=c(0,6),
ylim=c(0,6)
)
Overall, there seems to be an acceptable level of heterogeniety given the samples that are available, and there are no studies that contribute both to heterogeniety and influence on the overall result apart from Juvela.
Funnel plots - exploratory for small study publication bias
Studies with higher effect sizes are more likely to be published than those with lower effects. These missing studies which are not published are not identified and not integrated into the meta-analysis. This leads to publication bias meaning that the calculated effect size might be higher, and the true effect size lower since studies with lower effects were not published.
In addition, large studies are more likely to get published regardless of the effect size. Small studies are at the greatest risk, since they are only generally published if there is a large effect size.
Thus when assessing for publication bias, conventional assessment is focused on identifying whether small studies with small effect sizes are missing or not.
This can be performed using a funnel plot, however, we will later discuss why a funnel plot analysis may not be required.
funnel(pes.summary.glmm)
Here we can see that the funnel plot is assymetrical with visual assessment. The assymetry is primarily driven by 1 study in the top right corner with a large effect. We can identify the study by labelling the funnel plot.
funnel(pes.summary.glmm,
studlab = TRUE
)
There are many possible sources of funnel plot assymetry - selection bias from publication and selective outcome reporting bias, poor methodological quality leading to spuriously elevated effecs in small studies, true heterogeniety, statistical modeling artefacts and chance.
In our case, the most likely explanation for this appearance of the funnel plot is due to the outlier rupture risk in the Juvela study that contributes to true heterogeniety from true differences in rupture proportion.
In general, testing for funnel plot assymetry should always be performed in the context of visual assessment, and while there are many potential statistical tests for funnel plots including Egger 1997, Harbord 2006, Peters 2006 or Rücker 2008, even Cochrane suggests that tests for funnel plot asymmetry should be used in only a minority of meta-analyses (Ioannidis 2007b).
Sensitiity analysis
plot(inf.analysis, "es")
Here we rank the effect on the pooled rupture proportion using the leave-1-out method
Exploratory analysis
How do we assess the effect of exposure to prior SAH for participants that arise from the same population ?
Consider the data in the form of a 2 x 2 table, prior SAH as the exposure and rupture as the outcome.
ai = prior SAH and rupture +ve
bi = prior SAH and rupture -ve
ci = no prior SAH and rupture +ve
di = no prior SAH, and rupture -ve
n1i = ai + bi = total aneurysms with prior SAH
n2i = ci + di = total aneurysms without prior SAH
Rupture of the aneurysm is considered a rare event ie <1%, and the data are sparse with single 0s or double 0s in the 2 x 2 table.
This is methodologically challenging, and the choice of meta-analyis method is important. As we discussed The most common methods of MA is the inverse variance method, using the DerSimonian and Laird random effects model.
The DL method calculates an effect size separately for each study, with the standard error. The effect sizes are then synthesised across studies. However, when one of the cells has a 0 which is common with rare events, the inverse variance cannot be used because the variances become undefined.
There are 2 options for correction: use of a continuity correction ie adding a fixed value usually 0.5 or using calculating the risk difference. However using a continuity correction leads to excess bias in the effect, and can influence the result and conclusions (Stinjen 2010). Risk differences have poor statistical properties with too wide intervals when events are rare, and are also not recommended (Bradburn 2007)
There are also issues on how to handle double 0 studies, since these may also carry some meaningful data due to their sample size (Keus 2009).
In summary, there is debate on what is the best model to meta-analyse rare events. Use of Peto’s method to estimate the OR is often suggested for rare events, since this includes single zero studies without a continuity correction. Double zero studies are excluded.
However, to use Peto’s method, the following three conditions need to hold true ie a rare event <1%, the exposure / treatment groups are balanced, and the effects are not very large (Cochrane handbook). Unless all three are true, then Peto’s method also may give biased results. In our dataset, the groups with and without prior SAH are unbalanced ie >1:3, so Peto’s method is not appropriate.
Alternatively,the Mantel-haenszel without zero cell continuity correction can be used for unbalanced exposure / treatment groups.(Efthimiou 2018) This method provides a fixed effects meta-anlysis so is best used in the absence of heterogeniety and does exclude double zero studies. In our dataset, there is heterogeniety, and thus a random effects meta-analysis is preferred.
A generalised linear mixed method (GLMM) model can be used for odds-ratio meta-analysis for rare outcomes, specifically by utilising a hypergeometric-normal (HN) model for the meta-analysis of odds ratios (Stinjen 2010). This is appropriate since the event aneurysm rupture is not a true binomial distribution, but in fact a hypergeometric distribution.
The hypergeometric distribution is best explained by sampling coloured balls in an urn. Hypergeometric distribution is sampling without replacement compared to a binomial distribution where there is sampling with replacement, and the probability of success is required.
If the balls are of different weight or size, ie one has a greater chance of being chosen, this is a noncentral hypergeometric distribution. The non central hypergeometric distribution can be of the Wallenius type or Fisher type.
Wallenius type is the biased urn mode, where balls are taken out 1 by 1. Fisher type occurs when the outcome is known, and the number of balls in the urn and their colour need to be calculated. For large samples with a common outcome, the binomial distribution is a reasonable estimate. However in populations that are small, or outcomes that are rare, such as in our dataset where certain aneurysms have features that make them more prone to rupure ie heavier weighted ball, thus each outcome influences the probability of the next event. Thus the non central hyperegeometric distribution is required.
Recent developments in statistical packages, including in R, make these computationaly intensive methods practical and feasible. The HN model performs well, with minimal bias, and satsifactory coverage of the 95% CIs with rare events.
Mete-regression can also be included easily by extending the model to include a study level covariate.
Using metafor to asess the potential effect modification of prior SAH using a GLMM
Structure the data in this format
Author This signifies the column for the study label (i.e., the first author)
ai = prior SAH and rupture +ve
bi = prior SAH and rupture -ve
ci = no prior SAH and rupture +ve
di = no prior SAH, and rupture -ve
n1i = ai + bi = total aneurysms with prior SAH
n2i = ci + di = total aneurysms without prior SAH
Assumptions
For num_anr is not known, then num (of patients) is brought forward assuming 1 patient has 1 aneurysm.
For aneurysms in size cohort with prior SAH, this is estimated by using same proportion of patients with prior SAH in the total observed cohort, and applying this to the number of aneurysms in that specific size cohort.
Use the rma function from metafor to meta-analyse log odds ratio using conditional logistic regression model with random effects meta-analysis. The conditional generalized linear mixed-effects model with exact likelihood (model=“CM.EL”) avoids having to model study level variability by conditioning on the total numbers of cases/events in each study. For the odds ratio, this leads to a non-central hypergeometric distribution for the data within each study and the corresponding model is then a mixed-effects conditional logistic model.
Only studies that included patients both with and without history of prior SAH are included.
res <- rma.glmm(measure = "OR",
ai = ai,
bi = bi,
ci = ci,
di = di,
data = dat.psahpos,
slab = auth_year,
model = "CM.EL",
method = "ML"
)## Warning in rma.glmm(measure = "OR", ai = ai, bi = bi, ci = ci, di = di, :
## Studies with NAs omitted from model fitting.## Warning in rma.glmm(measure = "OR", ai = ai, bi = bi, ci = ci, di = di, : Some
## yi/vi values are NA.res##
## Random-Effects Model (k = 7; tau^2 estimator: ML)
## Model Type: Conditional Model with Exact Likelihood
##
## tau^2 (estimated amount of total heterogeneity): 2.1416 (SE = 2.6016)
## tau (square root of estimated tau^2 value): 1.4634
## I^2 (total heterogeneity / total variability): 56.2154%
## H^2 (total variability / sampling variability): 2.2839
##
## Tests for Heterogeneity:
## Wld(df = 6) = 3.2790, p-val = 0.7731
## LRT(df = 6) = 14.6179, p-val = 0.0234
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.1719 0.8681 0.1980 0.8430 -1.5296 1.8734
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1View the results as a forest
forest(res)
Overall, taking into account the limitations of the data, prior SAH may not increase rupture risk in small unruptured aneurysms measuring 10 mm or less. In addition, the overall rupture risk for conservatively managed aneurysms is more uncertain than previously considered with up to 2% risk of rupture.
These synthesised data analyses are considered exploratory and hypothesis generating, and additional data are required.
Exploratory - try and put the size subgroups all on one plot
additional code suggested by Thanh,
however, unlike retrospective / prospective subgroup analyisis, the subgroups are not mutually exclusive when conidering poportion of patients included with various size categories, ie numbers of 3 mm and less are included in 5 mm and less and so on. Thus I think that individual plots for size subgroup analysis submitted as supplementary files are more appropriate.
we can do additional subgroup analysis and combined plots for study quality etc.